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SYNOPSIS 


In India since the advent of river valley projects, 
irrigation through canal systems has brought about significant 
changes in the groundwater levels in the aquifers in canal 
command areas. These changes in hard-rock areas are so 
pronounced that drastic rise in water levels over the years 
has resulted in some localities through excess seepage from 
canals and less abstraction from wells. Thus, some of the 
regions in hard-rock terrains, once drought-prone, 


are 



currently facing or likely to face the threat of water-logging 
within the canal command areas. The consequences of the 
dominant seepage becomes apparent with time. The aquifer 
behaviour in such situations can be well understood, if all 
the inflows and outflows involved in a system are considered. 
The complex nature of the interaction between such features 
necessitates the use of mathematical modelling of aquifer 
behaviour. Once all the Inputs to the model are available, any 
of the parameters can be established using either the direct 
problem or the inverse problem. For example, the spatial and 
temporal distribution of water table heads can be worked out 
through direct problem using all other relevant input 
parameters. Similarly, in hard-rock areas, where information 
on aquifer parameters is scarce due to inadequate well testing 
facilities and procedures, estimation of these parameters can 
be done by modelling using inverse problem when all other 
system elements Including the spatial and temporal 
distribution of water table heads are available. If the 
information on all the system elements is known in a canal 
command area, it is not only possible to predict the water- 
logging at different places but also to suggest corrective 
measures through proper additional abstractions etc. The topic 
for the present study has thus been chosen keeping in view the 
problems of a canal command area. 

The study area is a part of the Krishna river basin and 
situated within the Nagarjuna Sagar left canal command area in 



Andhra Pradesh. Covering 830 sq. km, this area has the 
Nagarjuna Sagar left canal in the north and three rivers (Musi 
on west, Paleru on east, and Krishna on south) as its 


boundari es . 

Irrigation Is 

carrl ed 

for 

9 months in a year 

through a 

network of 

canals 

and 

abstraction through 


groundwater. Since the inception of the canals, regular 

monitoring is being done on the rainfall, water table heads, 

canal seepage and evaporation within the area and historical 

data are available for the same with the A.P. State 

Groundwater and Irrigation Departments. In the present work 

data pertaining to a period of ten years (1978 to 1988) at 69 

observation wells randomly distributed in the area were used 

for the modelling. The area has been divided into a 20 X 20 

2 

matrix with a grid size of 2. 5X2. 5 km and all the observation 
points are located on the same. The water table head data were 
processed to remove noise and the second spatial as well as 
the first temporal derivatives were calculated. These values 
together with the values for the well draft, Irrigation effect 
and rainfall recharge coefficient were assigned at each nodal 
point as inputs to a discretised Boussinesq’s equation and the 
aquifer parameters (transmissivity and storativity) are 
calculated for the entire period ( 1978 to 1988) in two blocks 
of five-year periods each with the help of a non-linear 
optimization routine. To evaluate the accuracy of the model, 
the water table head values were back-calculated using the 
estimated aquifer parameter values at each of the observation 
points. The deviation of the estimated head values from the 



actual head values was explained on the basis of the 
assumptions in the calculation of inputs as also the causative 
geological features. The accuracy of the model has also been 
tested by split sampling method. In addition, prediction of 
the changes in the water table situation at the end of next 
five year period (1993) was carried out using this model by 
incorporating a drainage channel and additional abstraction of 
ground water through wells. This pilot study was taken up as 
an illustration of the application of the model for water- 
logging problems. 

The work presented in the thesis is organized in six 
chapters. In Chapter 1, the topic has been introduced and 
details of the study area are given. The geological status of 
the area has also been presented. Unclassified granites of 
Archean age occurring in northern and central parts of the 
area are the dominant formations covering around 70* of the 
terrain. A few dolerite dykes with a general NE-SU trend are 
present In the western half of the granitic terrain. The 
granites are in turn overlain by quartz arenites and 
limestones, both belonging to Kurnool Super group (Upper 
Proterozoic in age). A few lineaments and two faults located 
on the basis of satellite imageries conform in their 
directions to the regional tectonic trend. The groundwater 
situation in the study area is described In this chapter and 
typical well hydrographs for all the three formations are 
presented. The scope of the present study has also been 



briefly dealt with. In Chapter 2, a review of modelling 
techniques for groundwater resource evaluation has been 
presented. The inverse and direct problems have been 
introduced. As the main emphasis was on the inverse problem in 
the present study, different methods of solving the inverse 
problem are presented with special reference to equation error 
criterion and the output error criterion. 

The processing of water table data forms the theme for 
Chapter 3. The relevant details of the methodology for water 
table head reconstruction using least square polynomial 
approximation are presented. The head data for sixty months 
(1978 to 1983) for all the 69 observation wells have been 
processed. The coefficients of the ploynomial fitted for this 
data have been calculated by considering an optimal form of a 
fourth degree polynomial equation which has been achieved by 
minimizing the standard error on considering a third degree 
polynomial. The decision for the same is made on the basis of 
students’ t test. The goodness of fit of the estimated 
coefficients has been checked through a multiple correlation 
coefficient. The estimated polynomial coefficients were 
utilized in computing the second spatial and first temporal 
derivatives, which were used later for the estimation of 
aquifer parameters. Typical water table contour maps and trend 
surface plots have also been presented. Effects of dykes and 
the geological structures (lineaments and faults) on the water 
table configuration have been discussed. 



In Chapter 4, the recharge and discharge components for 
the study area have been assigned for 120 months and at 69 
observation wells. A coefficient of 0.15 for recharge due to 
rainfall has been taken to be characteristic of the area. 
Baaed on the monthly evaporation rate obtained from A. P. State 
Ground Uater Department, losses due to evaporation were 
calculated. Recharge due to canal seepage has been estimated 
by considering conveyance losses for both the main and branch 
canals. 5% of the canal discharge has been considered to be 
Irrigation effect. 

Estimation of aquifer parameters has been carried out 
using the Inverse problem approach. The methodology along with 
the relevant mathematical background is presented in Chapter 
5. The Boussinesq’s equation was dlscretised with respect to 
system elements (transmissivity in X and Y directions, 
storatlvity, angle between the principal permeability 
directions and rainfall recharge coefficients). Taking the 
recharge and discharge parameters in a dlscretised form and 
derivatives of the water table head, an objective function has 
been formulated to estimate the system parameters. This 
function was solved by minimizing the sum of the squares of 
the residual errors with the help of a non-linear optimization 
routine (Sequential Unconstrained Flinimizat ion Technique). The 
optimal values of the system parameters for both the five-year 
blocks (1978-1983 and 1983 to 1988) were estimated when the 
error criterion is satisfied. The transmissivity estimates are 



2 2 
around 50 to 60 m /day with a variation upto 430 m /day in 

areas with dykes, lineaments and faults. The storativity 

estimates for the area are around 0.160, characteristic of 

unconfined aquifers. The angle between the principal 

permeability directions was estimated as 0.78 radians and the 

estimated values of recharge coefficient due to rainfall 

varied from 0.132 to 0.187. The groundwater storage has also 

been estimated for the area around each of the observation 

points . 

The accuracy of the model has been checked by simulating 
the head values, the estimated parameters and other related 
inputs using Alternating Direction Implicit Scheme. The 
deviation between the observed and estimated head values at 
different places was interpreted in terms of the geological 
factors such as the dykes and the assumption in the 
calculation of rainfall recharge. Effect of some of these 
factors has also been confirmed through use of split sampling 
procedure by taking estimated transmissivity values along with 
the recharge coefficients obtained for all the 69 observation 
points at the end of first five- year period as inputs for the 
estimation of heads for the second five-year period. In 
general, the error has been found to be less than 7% with 62 
out of 69 points having an error less than 5% . The values of 
transmissivity for all the 69 observation points for the two 
five-year periods along with the groundwater storage values 
have been presented in this chapter. 



The summary of the findings in terms of a comprehensive 
integrated picture is presented in Chapter 6. Application of 
the model for the study of water-logging problems has been 
illustrated. For this purpose, the groundwater configuration 
at the end of December 1987 has been chosen and water-logged 
areas have been Identified. Incorporating a drainage channel 
and additional drafts through wells, the aquifer behaviour has 
been simulated using the Alternating Direction Implicit 
Scheme. A reduction of water table up to 3 m In the 
water-logged areas has been observed in the predicted head 
distribution for the period ending May 1993. 



CHAPTER 1 


INTRODUCTION 


1.1 GENERAL 


Uith the Inception of Nagarjuna Sag ar project in Andhra 
Pradeah in 1966, considerable change haa cone about in the 
groundwater regine within the canal connand area. Canal 
irrigation haa reaulted in pronounced rise of the groundwater 
table at several placea in hard-rock terrain which haa 
hitherto been drought-prone with deep water table conditions. 
After two decades of irrigation practice, high water table haa 
created water-logging in several zones in the left and right 
canal connand areas. 

Prediction of aquifer behaviour as related to surface 
resources can effectively be carried out through nodelling. In 
the present work, a part of the Nagarjuna Sagar left canal 
connand area has been taken up for study to evaluate the 
groundwater situation. Data over a period of ten years (1978 
to 1988) pertaining to various inflows and outflows and water 
table heads were utilized to estlnate the aquifer paranetera , 
which in turn have been used to work out the groundwater heads 
at different periods of tine. The nodel has also been utilized 
for predicting changes in groundwater table heads by 



simulating additional draft through walls in arsas of high 
watsr table through an additional dralnag* channel 
incorporated in the model. 

1.2 THE STUDY AREA 

1.2.1 Location and Land Use 

The area chosen for study is situated in Nalgonda 
district of Andhra Pradesh and covers 830 sq.km with a maximum 
length and width of 45km and 60km respectively. Confined 
within latitudes 16°38N to 17°13N and longitudes 79°38E to 
80°08 E, the area has a general gradient of 1.75m/km, from 
north to south with highest and lowest elevations of 140m and 
60m respectively. Uith a good net work of fair weather roads 
providing accessibility to all villages, it is also well 
connected to the two prominent cities Vijayawada and 
Hyderabad. Of the total area of 830 sq.km of the basin, about 
500 sq.km is under cultivation. The rest of the area is 
generally barren and has scanty vegetation. The basin as a 
whole is deficient in rainfall. The area is irrigated with 
both canal water and groundwater for 9 months in a year. Rice 
cultivation is mainly restricted to areas adjacent to canals. 
In addition to rice, pulses and groundnut are the important 
crops raised (Table 1.1). On the whole. Irrigated area is 
about 60 percent of the total area of the basin. 

1.2.2 Hydrometeorology 

The basin has two rain gauge stations located in 
Huzurnagar and Kodad towns (Fig. 1.1). The data collected from 



Table 1.1 : Cropping Pattern in the Study Area 


crop 

area (' 

Rice 

51.00 

Groundnut 

15.20 

Chillies 

4.42 

Cotton 

3.53 

Millets 

0.80 

Vegetables 

10.70 

Pulses and Oil seeds 

13.35 


Sugar cane 


1.00 



thee* two places are utilised in computing the water budget. 
In addition, the evaporation data used in the study waa 
collected from the climatological station located at 
Huzurnagar. The annual evaporation is around 1800 mm. Based on 
the point rainfall measurements collected from the two rain 
gauge stations, annual rainfall distribution graph is prepared 
for the years 1978-1988 (Fig. 1.2). 

1.2.3 Surf ace Hydrology 

The study basin is bounded by three rivers and a canal. 
The Nagarjuna Sagar Left Canal (NSLC) which forms the northern 
boundary was commissioned in 1966, takes off from a dam on 
Krishna river at Nagarjuna Sagar with its flow from west to 
east. This canal with a length of 37km in the study area has 
several branch canals traversing the region. The hydraulic 
particulars of main canal are presented in Table 1.2. The area 
irrigated under each branch canal is shown in Table 1.3. 

The Paleru river on the eastern margin is controlled by 
canal supplies. The paleru lake at the extreme north of the 
basin is connected with a regulator to NSLC main canal and 
Paleru river. This lake acts as a balancing reservoir for the 
whole study area. The Nusi river on the western boundary has 
meagre flow and is filled with silt for maximum period in a 
year. This river has its confluence in the south with Krishna 
rlverwhlch has its flow from east to west and acts as a 
discharging boundary for the study region. 
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Pi^. i . 2 Rainfall distribution in the study area 






Table 1.2 : Hydraulic Particulars of Left Main Canal Uithin 

the Study Area 


sl.no. 

block 

no . 

1 ength 

Ckm) 

availability 

of 

water (cumecs) 

1 . 

8-11 

9.850 

225.17 

2. 

11-13(P) 

7.400 

169.46 

3. 

13(P)-15 

20.945 

157.60 


Table 1.3 : 

Hydraulic 

Particulars 

of Branch 

Canals 

al . name of 

block 

length 

discharge 

area irrigated 

no . canal 

no. 

(km) 

(cumecs) 

(»*) 


1. Nagulapadu 

8 

3.04 

0.3925 

35226.72 

2. Annaram-I 

9 

0.33 

0.1064 

817813.77 

3. Annaram-I I 

9 

1.25 

0.3579 

3186234.8 

4 . Somavaram 

9 

4.38 

1.0113 

9251012.2 

5. Chillapalll-I 

9 

9.47 

1.8501 

16943320 

6. Chi llapal 11- 1 1 

10 

7.55 

1.5322 

13805668 

7. Jonapadu 

10 

24.00 

15.7221 

134793000 

8. A.R.Gudem 

10 

9.11 

1.4791 

13526316 

9 . Ponugodu 

11 

2.27 

0.3534 

2834008 

10 .Sarvaram 

11 

8.74 

2.7012 

24375789 

11 .Kothaguda 

11 

1.34 

0.1192 

1093117.4 

12 .Kalavapalli 

12 

8.51 

2.0945 

18704453 

13 . J. P .Gudem 

12 

3.59 

0.6567 

5886639 

14 .Baithavole 

13 

3.62 

1.0294 

8910931 

15 .Barakatgudem 

13 

1.20 

0.3415 

1720648 

16 .13-D 

13 

1.50 

0.0872 

858299 

17 .Mutyala 

13 

29.41 

49.5829 

71457000 

18 .Komar abanda 

13 

8.61 

1.1475 

9943320 

19 .Polavaram 

14 

12.89 

2.3570 

3522267 

20 .Kotthagudem 

14 

26.81 

7.6563 

67797571 

21 .Raju Pet 

15 

3.89 

0.7717 

20109312 


* Source of data: A.P. State Irrigation Department. 
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1.3 GEOLOGICAL SET UP AMD GROUNDWATER SITUATION 
1.3.1 Geology 

A substantial portion (about 70 percent) of the study 
area has unclassified crystalline outcrops of Archean age 
composed mainly of granite with small pockets of hornblende- 
and quartz-schists. These are succeeded by the quartz arenites 
(of the Banganapalle group) which are in turn overlain by the 
Narji limestones (of the Jammalamadugu group), both belonging 
to the Kurnool Super Group (Fig. 1.1). The geological 
succession is as follows: 

Upper Kunool Super 

Group 

Proterozoic 

Unconformity 

Archean Unclassified Granites with basic 

Crystallines intruslves and schists 

The granites are basically of two types-the older grey 
type followed by the younger porphyritlc type. The grey series 
have lithological variation to pink granite. The porphyritlc 
division occurs as Intrusive laccollthic masses and dyke-like 
injections in the former type. The granites in the central 
region are highly weathered. Erosion along weak zones is quite 
prominent in the area south of Nagarjuna Sagar left canal and 
NE of Kamalacheruvu. Fine-grained dolerite dykes occur in 
granite and these are most Important in controlling 
groundwater condition of the study area. The strike of 
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dolerlte dyke a varies from NE-SU to NNE-SSU. These dykes 
follow prominent ve&k zone directions. The associated schists 
vary in composition from a quartz-schist to hornblende-schist. 
Ferrugeneous quartz-schists occurring as a patch SU of 
Betahvolu village have a trend N295° with foliation dipping 
towards NE. The quartz-hornblende schists located west of 
Lingagiri have a general strike of N320° with foliation 
dipping towards ENE. These schists are traversed by 
muscovite-bearing quartz veins. Details of lithology of 
formations encountered in a bore hole upto a depth of about 
70m in the Kodad area are indicated in Table 1.4. 

The sedimentary formations in the northeastern and 
southern parts of the area constitute a part of the Palanad 
sub-basin on the extreme northern part of the Cuddapah Basin. 
The quartz arenltes belong to the Banganapalle group and are 
coarse to medium grained with ferrugeneous matrix. The general 
strike of the formations is around N 80°-85° with low dips (5° 
to 8°) towards south. These formations exhibit current bedding 
of trough type with frequent changes in current directions as 
evidenced in some exposures. Indicative of their deposition in 
a shallow transgressive environment. The Narji limestone 
formations of Jammalamadugu stage occuring in the study area 
are located on the northern border of the Krishna river. These 
are represented by massive fine-grained limestones with 
occasional argillaceous pockets. The limestones are siliceous 
and have a NE-SU strike with varying dips towards south. At 
places the dips are very low in amount (about 5°) and the 
formations are horizontally bedded in the area NE of 



Table 1.4 : Details of Drill Hole Log at Kodad Town. 

(based on data form A . P . S . G . U . Department , A. P ) 


Formation 


Depth (m) Thickness(m) 


Top surface soil , lat erit e 

fine to medium-grained upto 3.0 3.0 

Granite, fresh, pink and grey with 
negligible weathering and fractures 

at 8.10 and 8.90 m 3.0 to 9.0 6.0 

Granite, pink with grey feldspar, 
weathered and fractured with clay 

fillings 9.0 to 11.0 2.0 

Granite, grey with pink feldspars and 
tnaflca, hard with negligible weathering, 

fractures at 18.8 and 20.0 m 11.0 to 24.0 13.0 

Granite, pink and grey, medium hard 

fresh but with fractures at 36.5 m 24.0 to 51.0 12.0 


Dolerlte, brittle and fine grained 51.0 to 70.0 


19.0 


Chltryala. 

The change in the course of the Krishna river at places 
in the study area appears to have been controlled by 
structural features. For example, the trend of the river in 
the northeastern part corresponds to the lineament (trending 
NNE-SSU) observed in the satellite imagery. Similarly the 
change in the river bend, after its confluence with Huai 
river, coincides with one of the two faults indicated in the 
geological map of the Cuddapah Basin (Geological Survey of 
India, 1981). These faults trending NE-SU extend into the 
study area (Fig. 1.1), one of them traversing the limestone 
while the other occurring in the quarts arenite formation. 
Uhile extensive fracturing and disturbance is evidenced within 
these formations in the vicinity of the faults, the 
dislocation is not seen in the field since these are 
intraf ormational faults and the outcrops are disturbed. The 
lineaments, as observed on the satellite imagery, have NE-SU 
trend. It is Interesting to note that the trend of these 
lineaments is similar in direction to that of the dykes 
present in the study area, possibly conforming to the regional 
tectonic trends. The tectonic framework of the region south of 
the study area is reported by Drury and Holt (1980). 

1.3.2 Groundwater Occurrence and Aquifer Characteristics 

Groundwater occurs under unconflned conditions. Ueathqred 
and fractured horizons constitute the aquifer zone. The wells 
in the area are both of large diameter type and tube wells. 
Their depths vary in different formations. In the granites and 



quartz arenites the wella are upto 13m while in the limestone 
area, the depth varies up to 22m. The water table head varies 
from around 144m in the north to around 54m in the south east 
corner of the study area. Typical hydrographs for for a period 
of 2 years (Jan. 84 to Dec. 85) for selected wells in all the 
three formations are presented in Figs. 1.3, 1.4 and 1.5. 
Transmissivity (T) and storativity (S) values for the granites 
as reported by the Andhra Pradesh State Groundwater department 
are in the range of 60-170m*/day and 0.2-0.24 respectively. No 
information is available on the aquifer parameters of the 
sedimentary formations in the study region. 


1.4 SCOPE OF THE PRESENT STUDY AND ITS ORGANIZATION 

The groundwater regime in an area is controlled by 
several parameters such as recharge from rainfall and canal 
seepages, aquifer characteristics, groundwater draft and flows 
across the boundary. A proper understanding of the interplay 
of all these parameters can be achieved through modelling. 
Modelling forms an integral part of groundwater 
investigations. In the present work a portion of left canal 
command area of Nagarjuna Sagar project in Andhra Pradesh has 
been chosen for the study. Surrounded by three rivers and a 
canal as its boundaries on four sides, the area has been 
irrigated through canal irrigation for over two decades. 

The data for a period of ten years between June 1978 to 
May 1988 pertaining to various components Involved have been 
utilized in the model. Processed water table head values after 
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removing the noise in the data have bean used together with 
the other calculated recharge and discharge values as inputs 
and the aquifer parameters (transmissivity and storativity) 
have been estimated using the inverse method. These estimates 
were in turn utilized for checking the validity of the model 
through direct problem approach by simulating the water table 
head data, which were compared with observed field data. The 
model has also been applied to monitor the water table changes 
with additional well drafts Incorporated. 

The material contained in the thesis is organized and 
presented under several chapters. 

In Chapter 1 the study area is Introduced in terms of its 
location, land use, meteorological background, geology and 
groundwater occurrence. The scope of the present study is also 
outlined. 

In Chapter 2, a review of modelling techniques for 
groundwater resource evaluation has been presented. The 
inverse and direct problems have been Introduced and their 
scope indicated. 

The processing of water table data forms the theme of 
Chapter 3. The relevant mathematical background for water 
table head reconstruction using least square polynomial 
approximation has been presented in detail. The processing of 
water table data for the study area has been presented. Two- 
dimensional (contour maps) and three-dimensional (trend 
surface maps ) graphical representation of the water table 
head variation for different periods of time are presented. 



Estimation of rtch&rge and discharge components and the 
recharge evaluation for the rainfall data together with the 
estimation of groundwater system elements are developed in 
Chapter 4. The recharge and discharge components for the study 
area are estimated for 120 months and at 69 observation 
points . 

Chapter 5 deals with estimation of aquifer parameters 
using the inverse problem approach. The methodology with 
relevant mathematical background pertaining to discretisation 
of governing equation, estimation of derivatives and 
optimization has been presented. The aquifer parameters and 
the groundwater storage have been worked out for the study 
area. In addition, the rainfall recharge coefficients and the 
angle between principal permeability directions were obtained. 
The calculated aquifer parameters have been used for testing 
the validity of the model. Estimation of water table heads has 
been carried out by simulation and the comparison of the 
estimated values with the observed data has been attempted to 
analyze the error range. The anomalous values of the water 
table heads have been interpreted in terms of basin 
charact eristics . 

An integrated picture has been presented in Chapter 6. 
The application of the model for monitoring the groundwater 
variation with additional well draft in areas of shallow water 
table has been also illustrated by simulating the situation 
for a particular set of values. 



CHAPTER 2 


GROUNDWATER RESOURCE EVALUATION BY MATHEMATICAL MODELLING 

- A REVIEW 


2.1 INTRODUCTION 

Significant research work has been conducted in the field 
of groundwater modelling in the last three decades. It is used 
extensively in understanding aquifer flow, as well as in 
evaluating regional groundwater resources. If adequate data is 
available, it is an ideal tool for development of basin-vise 
plan. A review is attempted in the present chapter to bring 
out the main applications with emphasis on inverse problem 
pertinent to the present study. 

2.2 MODELLI NG TECHNIQUES 

Model is a tool designed to represent a simplified 
version of reality (Uang and Anderson, 1982). A mathematical 
model consists of a set of linear differential equations that 
are known to govern the flow of groundwater. Mathematical 
models of ground water flow have been in use since the 
beginning of this century. The reliability of predictions 
using groundwater model depends on how well the model 
approximates the field situations. Since the field situations 



are too complicated to be simulated exactly, assumptions must 
always be made In order to construct a model. Usually such 
assumptions necessary to solve a model analytically are 
restrictive . 

One can consider two types of models-f inite difference 
models and finite element models. In either case, a system of 
nodal points are superimposed over the problem domain. In the 
finite difference models, nodes are located either at the 
centre of the grid lines as in a rectangular grid pattern 
(Fig. 2. la) or inside the cell as in polygonal grid pattern 
(Fig. 2. lb). Uhile in the former case mesh-centered nodes are 
used, block- centered nodes are in use in the latter. The 
concept of elements is fundamental to the development of 
equations in the finite element method. 

In finite difference modelling, the governing equation of 
flow is approximated using finite differences and resulting 
set of linear or nonlinear algebraic equations is solved by 
using direct or iterative techniques. Stallman (1956) was the 
first to apply numerical methods to groundwater hydrology. 
Aquifer problems in which transmissivity is independent of 
variations of the piezometric head are linear problems. Thus 
confined and leaky aquifers generally lead to linear problems 
while unconfined aquifer problems may be nonlinear. 

For unconfined aquifers, the solutions are greatly 
facilitated if Dupit-Forchheimer "a assumptions are considered. 
In such a case, the governing equation (Remson et al ,1971) for 
two-dimensional transient groundwater flow in an isotropic, 
homogeneous, unconfined aquifer, with no vertical accretion is 




(a) Rectangular grid 



(b) Polygonal grid 

FIG. 21 GRID PATTERNS FOR FINITE DIFFERENCE 
MODEL 
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where K is the permeability, h is the hydraulic head, x 
and y represent distances in the respective directions of 
flow, t is the time increment and S is the specific yield. 
Taking into consideration the anisotropy, nonhomogeneity , and 
the vertical accretion, Eqn.2.1 has been modified (Hollond and 
Dracup, 1970) as 
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where x and y refer to the principal permeability 
directions, and and K yy are the permeablltles in 

directions x and y respectively and Q is the net inflows and 
outflows with change in storage. 

The governing differential equation (Rushton and Rathod, 
1985) for three-dimensional time-variant flow in an aquifer is 



(2. 3. a) 

where K , K and K are the hydraulic conductivities in 
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In which z is the height above an arbitrary datura, p(x,y,z,t) 
is the pressure, p is the density of the fluid and g is the 
acceleration due to gravity. 

2.3 AQUIFER RESPONSE 

The water balance equation proposed by Sokolov and 
Chapnan (1974) can be written as 

Inflows - Outflows - Increase in storage + Error tern = 0 

...(2.4) 

« 

Inflows Include the rainfall recharge, recharges from 
rivers, subsurface horizontal inflow, artificial recharge and 
inflows frora other aquifers (overlying or underlying). Out- 
flows include baseflow to the rivers, outflow of groundwater 
into the zone of aeration for noisture recovery lost by 
evapotranspiration, outflow to overlying and underlying 
aquifers, subsurface horizontal flow, groundwater flow through 
springs and groundwater punpage. The response of any aquifer 
in terras of Inflows and outflows can be evaluated through the 
lunped models (Sokolov and Chapnan , 1 9 7 4 ) or distributed 
nodela (Reason et al, 1974; Plnder and Gray, 1977). The 
lunped aquifer response to known Inflows and out- flows can be 
obtained by designing groundwater storage fluctuations in 
terras of groundwater head fluctuations and storatlvity. 



The distributed groundwater flow models are based upon 
the solutions of differential equations governing two- 
dimensional (Eqn.2.2) or a three-dimensional (Eqn. 2.3a) case 
like the representing transient groundwater flows in saturated 
zone. Closed form or series solutions of governing 
differential equations are available for only Idealized 
boundary and recharge conditions (Be*?, 1967). These are 
generally based upon the assumptions of homogeneity and 
isotropy. For predicting the aquifer response to the spatially 
and temporally varying recharge and pumpage under realistic 
conditions, the governing differential equations have to be 
solved by appropriate computer-assisted numerical methods. The 
discretisation of space, necessary for finite difference 
approximations, can be based upon either any convenient 
pattern (Tyson and Uebor, 1964) or a more stringent 
rectangular grid pattern. The former is generally known as 
'Tyson and Uebor' model and the latter as 'finite difference 
method' . 

The finite difference approach Introduced by Richardson 
(1910) Involves calculation of approximate solution of partial 
differential equations. This method has been programmed to 
solve two-dimensional (Bittenger et al , 1967; Finder and 
Bredehoeft, 1968; Pinder, 1970, Prickett and Lonnquist, 1971; 
Bedinger et al , 1973, Trescott et al, 1976) and three- 
dimensional (Bredehoeft and Pinder, 1970; Trescott, 1975) 
transient groundwater flow problems. The implicit form of the 
finite difference equation generally requires the solutions of 
large number of linear simultaneous equations and hence the 



requirement for storing all the coefficient matrices is 
enormous. This difficulty is overcome by employing various 
iterative procedures which primarily bank upon the sparseness 
of the coefficient matrix. 

Some of the more commonly used procedures are alternating 
direction implicit scheme (Peaceman and Rachford, 1955), 
Crank- Nlcolson scheme (1947), line successive over relaxation 
(Trescott et al, 1976) and strongly Implicit procedure (Stone, 
1968). The numerical properties of these methods have been 
extensively studied (Rushton, 1974; Tomllson and Rushton, 
1975; Trescott and Larson, 1977). Murray and Johnson (1977) 
have demonstrated the second method of llnerlsation of 
Bousalneaq’s equation for simulating the response of 
unconfined aquifers. The finite element method is reported to 
overcome many difficulties relating to irregular geometry of 
the area, heterogeneity and boundary conditions (France, 
1974), in addition to giving results of higher degree of 
accuracy (Javandal and Uitherapoon, 1968). Other significant 
contributions are by Neuman and Uitherapoon (1970), Pinder and 
Frind (1972), Neuman (1973a), Pinder (1973), Pinder and Gray 
(1977). 

2.3.1 Recharge 

Recharge rate, an input requirement of all aquifer 
response models, is one of the Important factors controlling 
the amount of water that may be pumped from an aquifer without 
depleting it. Rushton (1986) has considered the following 
equation for describing flow in an aquifer. 
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where h is the groundwater head, 

T x and Ty are the transmissivities, 
x and y are the space coordinates, 
t is the time coordinate, 

S is the storage coefficient or specific yield, 
and Q is the recharge intensity 

The Richard’s equation for estimating groundwater recharge 
using finite difference method (Krishnamurthi , 1977; Kafri and 
Asher, 1978) is as follows: 
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c 

where — is capillary pressure, & is moisture content, z and 
t are spatial and temporal coordinates while K(A) is the 
hydraulic conductivity. 

The Input data requirement for groundwater modelling 
would amongst others include transient moisture storage 
measured as a function of vertical position, permeability-soil 
moisture curves and evapotranspiration. The approach for the 
recharge estimation is based upon water budgeting (Eqn.2.7) 
and field capacity. 


P = ET + SRF + ASM 


(2.7) 



where P is the precipitation, BT la actual evapotran- 
apiratlon, SRF la the aurface run-off and ASM la the increaae 
of water availability in root zone. Recharge (Eqn. 2 .8) la 
related to current aoll moisture storage at field level. 

Recharge - A SM - DFCT if DFCT < A SM 

= 0 if DFCT &A SM (2.8) 

where DFCT la the current aoll moisture deficiency below field 
capacity. 

Some doubts have been expressed regarding the validity of 
this conventional method (Uatson et al , 1976; Ruahton and 
Uard, 1979). Ruahton and Uard (1979) have proposed alternative 
recharge mechanism which allows recharge to occur even when a 
aoll moisture deficit exiata. Morel-Seytoux (1976) presented 
formulae for prediction of pounding time and cumulative 
infiltration, under influence of rainfall. Reeves (1975) haa 
teated a method of recharge estimation, which assumes that the 
maximum Infiltration rate la simply a function of cumulative 
rainfall regardless of rainfall versus time history. A few 
reported studies make use of water table data for the 
estimation of recharge (Venetis, 1971). 

The experimental methods for the measurement of recharge 
Include tritium tracer technique (Vogel et al, 1974; Sukhrla 
et al, 1976) and gamma ray transmission (Singh and Chandra, 
1977). An empirical relation between the rainfall and recharge 
was given by Chaturvedl and Chandra ( 1961 ) . 



2.3.2 Interpolation of Groundwater Head 

The need of interpolation arises because the location of 
the points of groundwater data is very rarely decided on the 
basis of the requirements of the digital models. Thus, in most 
of the cases the groundwater heads at the nodal points have to 
be estimated by Interpolating the recorded heads. The 
conventional graphical and lagranglan methods of Interpolation 
(Ralston, 1965) may not be always suitable because of the 
subjective bias and artificial undulations associated with the 
high degree polynomials respectively. Sagar et al (1973, 1975) 
demonstrated the use of spline functions to approximate the 
spatial variation of groundwater heads. The closed form 
functions so obtained can assist in Interpolation. Delhomme 
(1978) has reported the possible use of kriglng technique in 
providing spatial variation of groundwater head. Blrtles and 
Morel (1979) proposed a least square based method of 
Interpolation wherein the unsteady state flow equation is 
reduced to a steady state by assuming zone-wise spatial 
uniformity of certain equations. 

Kashyap and Chandra (1982) have evolved the least square 
polynomials of spatial coordinates p and q in any two 
arbitrarily selected directions to approximate spatial 
variation of water table data elevations. 

2.4 L INVERSE PROBLEM 


The major hurdle in mathematical modelling of groundwater 
flow is the collection of adequate data relating to the 
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spatial and temporal distribution of transmissivity and 
storativity. Of the difficulties that arise, the principal one 
is that these properties are not physically measurable 
quantities (Yeh, 1975). These are, in fact, parameters 
appearing in the differential equation governing the aquifer 
response (groundwater head fluctuations) to the aquifer 
excitation (pumpage, recharge or change in boundary 
conditions). The estimation of transmissivity and storativity 
is an inverse process involving the historical excitation, 
response and associated initial and boundary conditions. The 
system in this case is represented as in Fig. 2. 2. A variety 
of^nethods exists based on the parameters employed for the 
estimation of the aquifer response to different forms of 
excitations with varying conditions and assumptions. 

The most widely used method is generally test pumping 
(Kruseman and Ridder, 1970). This method essentially consists 
of generating and recording aquifer response under a 
controlled excitation in the form of pumping from a single 
well. The generated data are employed to estimate parameters 
based upon the criterion of 'closest match’ between the 
recorded response and the response given by the governing 
differential equation. An accurate estimation of the 
transmissivity and the storativity through conventional pump 
test methods requires all the real field situations being 
properly incorporated in the equation involved. For example, 
Thels (1935) equation assumes axlsymmetric radial flow towards 
a fully penetrating discharging well of infinitesimal diameter 
in an infinite confined aquifer. Although, several modlfi- 




FIG. 2-2 DIRECT AND INVERSE PROBLEMS IN GROUNDWATER 
HYDROLOGY 
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cations have bean reported (Jacob, 1950; Plnder and 
Bredehoeft, 1968; Ualton, 1970; Yeh and Tauxe, 1971; Marino 
and Yeh, 1973;' Lakshminarayana and Rajagopalan, 1978), there 
la no single comprehensive method incorporating all the field 
condi tiona . 

Estimation of transmissivity and storativity can also be 
carried out using the fluctuations of water levels in 
hydraulically connected water bodies, along with other form of 
aquifer exication. These can be viewed as excitations induced 
by change in boundary conditions. Pinder and others (1969) 
employed iterative procedures to determine aquifer dlffualvlty 


from the 

historical 

data 

relating 

to the 

water 

table 

elevations 

and river 

stage. Yeh 

(1975) 

used 

quasi- 

linearization (Bellman 

and 
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diffuslvity of an unconfined aquifer baaed on historical data 
of stream-aquifer interaction. Singh and Sagar (1977) 
developed an analytical method to determine aquifer 
diffuslvity using the Green's function of linearised 
Bousalnesq's equation and specifying an extra boundary 
condition at the stream-aqulf er Interface. Sagar and Singh 
(1978) studied the effect of observational errors in raw data 
on the values of the aquifer diffuslvity, the most restrictive 
among them being one-dimensional flow, absence of vertical 
accretion and existence of Improper boundary conditions at the 
other end of the aquifer. These two methods can thus be 
employed to estimate aquifer parameters of flood plains. The 
study reported by Venetls (1971) overcomes the assumption of 
‘no vertical accretion’ in a limited way by using uniform 



accretion. Other reported contrlbutiona are by Ferris (1952) 
and Rove (1960). 

The estimation of aquifer parameters Involves the 
excitations, response data such as the fluctuations in 
groundwater heads in space and time in response to known 
pumping and/or recharge pattern, the initial and boundary 
conditions. Since the closed form analytical solutions for 
such complex system are not available, such an approach has 
essentially to be based on the Bousslnesq’s equation and is 
treated as 'inverse problem’ in groundwater hydrology. There 
are two types of errors associated with the inverse problem. 

1) The system modelling error, as represented by 
performance criterion, 

2) The error associated with parameter uncertainity 

(Teh, 1986). 

An increase in number of unknown aquifer parameters will 
generally Improve the system modelling error, but will 
increase the aquifer parameter uncertainity and vice versa. 
The optimum level of parameterization depends on the quantity 
and quality of data. The inverse problem is often ill-posed, 
due to non-uniqueness and instability. Small errors in the 
head will cause serious errors in the identified aquifer 
parameter leading to instability. 

2.4.1 Classification of Parameter Identification Methods 

Various techniques have been developed to solve the 
inverse problem of parameter identification. Neuman (1973) 



classified the techniques into two categories- 'direct * and 
'indirect 1 . The direct approach treats the model parameters as 
dependent variables in formal inverse boundary value problem. 
The indirect approach is based upon the output error criterion 
where repetitive numerical solution of Boussinesq’s equation 
with successively modified aquifer parameter values are 
calculated to reach the closest reproduction of recorded 
aquifer response under the recorded excitations (withdrawals 
and recharges) and initial as well as the boundary conditions. 
The modifications in the aquifer parameter values at the end 
of each iteration can be baaed either upon the mathematically 
rigorous procedures (Knowles et al , 1972) or upon subjective 
judgment. The latter approach is also known as 'calibration* 
or 'subjective optimization*. 

Kubrusly (1977) classified parameter identification 
procedures of distributed parameter system into three 
categories : 

1) The direct method, which consists of those methods 
that use optimization techniques directly to the distributed 
parameter model 

2) Reduction to lumped parameter system, which consists 
of those methods that reduce the distributed parameter system 
to a continuous or a discrete time lumped parameter system 
described by ordinary differential equation or difference 
equation 

3) Reduction to algebraic equation, which consists of 
those methods that reduce the partial differential equation to 
an algebraic equation. 





2.5 DIRECT METHOD USING THE EQUATION ERROR CRITERION 

In the direct methods of solving Inverse problem, 
Boueslnesq's equation Is directly solved for the aquifer 
parameters. In practice, observation wells are sparsely 
distributed In the flow region In an arbitrary fashion and 
only a United number of observation wells are available. To 
formulate the Inverse problem by the equation error criterion, 
missing data have to be estimated by Interpolation. The 
Interpolated data along with the observations, which also 
contain noise, on substitution In the governing differential 
equation results In an error term (known as the equation 
error). This error Is minimised over the proper choice of 
parameters (Teh, 1986). 

The equation (Jacob, 1950) for confined aquifer which 
Incorporates anisotropy, heterogeneity and vertical accretion 
Is expressed as 
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where and T are transmissivity values In the 
direction of principal permeabilities collinear with x and y 
and 0 is the net vertical accretion rate per unit area. The 
above equation, rewritten with transmissivity as unknown, 
becomes a first order partial differential equation and can be 
expressed In the form 
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To estimate transmissivity (T) values from this partial 
differential equation, values of T at the boundaries 
designated as Cauchy data (Emaellem and delfarslly, 1971) are 
needed. This approach can be employed in the estimation of 
spatially distributed transmissivity values from the 
groundwater data and vertical accretion data of single time 
period, provided the first derivatives of groundwater data are 
elaborate enough. Since a single uncoupled equation can yield 
the solution of one variable only, this method can not be used 
to calculate both transmissivity and storatlvlty (S). Either 
the S or the groundwater data corresponding to steady state 
conditions has to be an input data. Since in most of the real 
situations the necessary Cauchy data are not available, this 
approach has its limitations. This problem can be overcome by 
employing multiple period data. The overdetermined system, so 
obtained can be utilized to arrive at optimum values of the 
parameters, by minimizing the predecided function of the 
residues in the Boussinesq's equation. This approach 
introduces the 'smoothening ' effects if any discrepancies are 
present in the raw data. 

2.5.1 Solution Procedures 

Various methods are available to solve inverse problem by 
equation error criterion (the direct method) such as energy 



dissipation method (Nelson, 1968), linear programming approach 
(Klelnecke, 1971), the use of a flatness criterion (Emsellem 
and deNarsily, 1971), the multiple objective decision process 
(Neuman, 1973), the Galerkln method (Frlnd and Plnder, 1973), 
the algebraic approach (Sagar et al , 1975), the Inductive 
method (Nutbrown , 1975), linear programming and quadratic 
programming (Hafez, 1975), minimization of quadratic objective 
function with penalty function (Navarro, 1977), and the matrix 
Inversion method allied with kriglng (Yeh et al , 1983). To 
minimize the instability and nonuniqueness, regularity 
conditions are often required. 

Sagar and co-workers (1973,1975) proposed the 'algebraic 
equation method* which dispenses with the requirement of 
boundary conditions and it Is not based on the discretisation 
of space by nodal points. The spatial derivatives of hydraulic 
head are estimated by one-dimensional spline function fitted 
to the observed hydraulic head data. The boundary conditions 
are eliminated by studying water budget at each individual 
space point and the parameters estimated by minimizing the sum 
of squares of residual errors. However this method Is useful 
for the estimation of transmissivity only. 

Emsellem and deNarsily (1971) for the first time 
demonstrated the use of the multiple objectives, in context of 
aquifer parameter estimation. One of the objectives Is to 
minimize the residue functional and the other is to affect 
'smoothness* with respect to variation of transmissivity In 
space. The suggested approach is to significantly reduce the 
residue functional gradually. The method applies to steady 
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state conditions and, thus, can yield only the estimates of 
transmissivity. The idea of multiple objectives has been 
further extended by Neuman (1973b, 1975), who has proposed a 
number of criteria to affect 'smoothness' solutions. The 
problem has been solved by applying linear programming to 
finite element-based model of groundwater flow. 

Nutbrown (1974) has proposed another viable method of 
estimating transmissivity and storativity. The estimation of 
transmissivity is done by the water balance of stream-tube for 
a period of zero storage change. The stream-tube is extended 
to a turning point (water divide) and the transmissivity at 
the other end is determined. The storativity is obtained by 
employing the data of the period by displaying the storage 
changes. This approach was further extended by Brltlea and 
Morel (1979) who incorporated 'equivalent steady state* 
concept for interpolation. Other significant contributions are 
the study of input of additional data on calibration by Gates 
and Kislel (1974), regression-based drawdown model by Haddock 
III (1976) and the subjective information by Lovel and 
co-workers (1972). Aquifer excitation along with 
transmissivity and storativity was estimated by Tlaon (1965), 
Monech and Kisiel (1970), Venetis (1971), and Sagar et al 
(1975). 

2.6 INDIRECT METHOD USING THE OUTPUT ERROR CRITERION 

The criterion used in this approach is generally the 
minimization of a 'norm* of the difference between observed 



and calculated heads at specified observation points (Yeh, 
1986). It has both advantages and disadvantages. While it can 
be applied with limited data and does not require 
differentiation of the measured data, the disadvantage of this 
approach la that the minimization is usually nonlinear and 
non-convex. The algorithm starts from a set of initial 
estimates of the parameters which Improves with each iteration 
until the system model response is sufficiently close to the 
observed data. 

2.6.1 Solution Procedures 

Gauss-Newton procedure for inverse problem solution was 
first applied by Jacquard and Jain (1965), followed by Jahns 
(1966). Along with quadratic interpolation, this method was 
used by Thomas et al (1972) and also adopted by Yoon and Yeh 
(1976a, 1976b), Yeh and Yoon (1981) and Sadeghipour and Yeh 
(1984) along with Rosen’s gradient projection. Conjugate 
gradient oriented method proposed by Gauss-Newton and further 
modified by Narquardt has been used by Gavalas et al (1976). 
Cooley (197.7) modified Gauss-Newton method using nonlinear 
regression by linearisation. Cooley (1982) applied this method 
with the help of prior information to solve Inverse problem. 
Shah et al , (1978) and Sun and Yeh (1985) have also used 
Gauss-Newton method to estimate aquifer parameters. 

Steepest descent method was first adopted for Inverse 
problem by Vemurl and Karplus (1969). Chen et al (1974) used 
the same method along with conjugate gradient. Chavent and 



co-workera (1975) have applied this method without conjugate 
gradient for oil reservoir. Conjugate gradient vaa alone 
applied to solve inverse problem by Neuman (1980). 

Quasi linearisation was first used by Yeh and Tauxe 
(1971) to identify aquifer parameters. Later Narlno and Yeh 
(1983) and Distefeno and Rath (1975) have applied the same 
method to solve Inverse problem. The other methods of solving 
Inverse problem worth mentioning are Newton-Raphson method 
used by Neuman and Yacowitz (1979), maximum likelihood and 
krlglng method of solution by Kitanldis and Vomvoris (1983), 
and cokriging by Hoeksema and Kitanldis (1984), Gaussian 
conditional mean and krlglng as compared by Hoeksema and 
Kitanldis (1985). Quasi linearisation, maximum principle 
gradient, influence coefficient and linear programming 
algorithms were compared by Yeh (1975a). Yeh (1975b) used 
quadratic programming to solve inverse problem. 

It can be proved that least square solution approach 
yields the most likelihood solutions when the errors in raw 
data are normally distributed (Beck and Arnold, 1977). This 
least square form of the residue functional, necessitates the 
use of nonlinear programming. Cooley (1977) proposed the 
method of regression for the parameters of two-dimensional 
steady state flow domain and also analyzed the goodness of 
fit and significance of computed parameters in a subsequent 
paper (Cooley, 1979). Kleinecke (1971) adopted 'mini-max' and 
the 'sum of the moduli’ functionals and used the linear 



programming to obtain optimal values of the aquifer 
parameters. The main drawback of this method is that the final 
solutions may not include the parameters of all the grid 
points. For example, out of 144 nodes, 44 did not generate an 
estimate of storativity. Results of transmissivity are even 
more sparse, since out of the 676 estimates, 56 percent were 
not in the solution domain. Subsequently , Navarro (1977) 
modified the Klelnecke’s method to make use of field 
information relating to the values of transmissivity and 
storativity. 

Haimes and others (1968) employed decomposition and 
multilevel optimization, to arrive at optimal parameters. This 
method assumes an infinite aquifer. Labadie (1975) proposed 
the term 'surrogate parameters' meaning that the estimated 
parameters may be quite different from the aquifer properties, 
since those parameters represent a sort of lumped effect of 
aquifer properties and many other unaccounted local features 
like heterogeneity and complex boundary conditions. 

Kashyap and Chandra (1982) have developed a numerical 
scheme to estimate quantitatively the parameters related to 
geohydrological and hydrological characteristics of aquifers, 
employing historic data like hydraulic head, rainfall and 
pumpage. This scheme is based upon constrained minimization of 
the sum of squares of residuals in the Boussinesq’s equation. 
Derivatives of the hydraulic head are estimated by polynomial 
approximation of least squares. 



Thus, the inverse problem for the estimation of aquifer 
parameters from aquifer excitation and other related initial 
and boundary conditions can be solved by using equation error 
criterion (direct method) or output error criterion (indirect 
method). In the present work, the direct method has been 
adopted for the estimation of aquifer parameters. As already 
indicated, groundwater head values at various nodal points as 
obtained by interpolation from available field data are needed 
as inputs. The procedure for the interpolation and smoothening 
of water table head data is outlined in the next chapter. 



CHAPTER 3 


WATER TABLE HEAD RECONSTRUCTION 

3.1 INTRODUCTION 

For quantitative evaluation of groundwater potential, 
historic record of water level data is a prerequisite. The 
available records usually consist of groundwater head measured 
periodically in observation wells (open wells or piezometers). 
These discrete point data can be employed to generate 
continuous functions which approximate the true functional 
relationships between the groundwater head and the spatial 
coordinates in the domain bounded by the problem area. Such 
approximate functional relationships can be greatly useful in 
subsequent analysis of the data, to estimate aquifer 
parameters. A deterministic method for generating such 
approximation has been developed in this chapter. The method 
consists of regression of groundwater head data by minimizing 
sum of squares of residual errors, with relevant statistical 
tests . 

< 

3.2 PROCESSING OF WATER TABLE DATA 

The observation well data can be processed to yield a 
variety of Information relating to the groundwater movement. 



storages and groundwater head elevations at various locations 
in the problem domain. There are three basic numerical 
techniques for this data processing viz., differentiation, 
integration and interpolation. 

The differentiation of groundwater head data with 
respect to space is necessary because of the linear velocity- 
hydraulic gradient relation incorporated in Darcy's law. Thus, 
the first derivatives of groundwater head with respect to 
space i.e. the hydraulic gradient needs to be calculated for 
the estimation of subsurface velocities and recharge across 
the boundary. Similarly, the second spatial derivatives of 
groundwater head assist in calculation of the imbalance 
between horizontal Inflows and outflows at any space point 
across the boundary. 

The integration of groundwater head in space is required 
to be carried out for estimating the water released from or 
taken into the storage in a given area for a certain period. 
The interpolation becomes necessary, when the nodal points of 
a distributed parameter groundwater flow model do not coincide 
with the observation points. The available data of the 
observation points are employed to assign groundwater heads at 
the nodal points. 

3.3 EXISTING METHODS 

The mathematical operations mentioned above can be 
performed by arriving at a spatial variation of groundwater 
head using stochastic, deterministic or graphical method. 
Krlglng, a stochastic method (Delhomme, 1978) assumes that the 



groundwater system Is too complex to be explained by 
analytical expressions. For deterministic system analysis, the 
spatial variation can be arrived at by graphical procedures or 
functional approximation. Alternatively, finite differences 
can be employed to perform the Interpolation, differentiation 
or integration. 

The graphical procedures essentially Involve the drawing 
of contours of equal groundwater heads by visual interpolation 
of the groundwater head elevations at the observation points. 
This procedure, though employed quite frequently, has two 
major handicaps in that a strong subjective bias and a 
subconscious assumption of linear variation of groundwater 
head between two observation points are Involved. These 
handicaps become all the more pronounced when the data points 
are sparse or non-uniformly distributed in space. The finite 
difference methods of differentiation. Integration and 
interpolation are procedures which, in most of the cases, 
afford the estimations of error bounds. However, the 
requirement of data along the two orthogonal directions is 
rarely satisfied. 

The 'functional approximation’ approach consists of 
approximating the true functional relation between the 
groundwater head and space coordinates by an explicit function 
which is directly amenable to the mathematical operations. 
Functional approximation can be either 'exact type’ l.e. no 
residues at the observation points or the 'least-square’ type 
(Ralston, 1965) with residues at the observation points. These 
procedures require the total number of coefficients in the 



approximating function to be equal to or less than the total 
number of data points. Uhere the general form of the 
functional relation is not known, a polynomial may be used to 
approximate the function in a domain, provided it is known a 
priori that the function is continuous. This approach is 
validated by the Ulerstr&ss’s theorem (Ralston, 1965). The 
lagranglan methods of "exact' functional approximation in one- 
dimension by a polynomial can be suitably extended to a 
two-dimensional situation. This method, however, requires a 
polynomial of very high degree which apart from increasing the 
computational efforts involved in the numerical operation of 
the function, can also cause large artificial undulations of 
the functional estimate in between the observation points. 

The spline functional approximation (Sagar et al , 1973 , 
1975) involves passing piece-wise continuous polynomial 
function through known functional values at the observation 
points with the compatlability conditions at the interfaces of 
the adjacent polynomials satisfied. This approach is one of 
numerical curve-fitting methods. It overcomes the necessity of 
employing higher degree of polynomial but may still cause 
artificial undulations if the observation well data are not 
completely error-free, as they would rarely be. The 
differentiation of these approximating functions to arrive at 
the derivatives of the true function, may involve large errors 
originating from even small amplitude of "noise' present in 
the raw groundwater head data. 

The magnitudes of noise are not known but the frequency 
distribution will be near Gaussian if the sample size is large 



enough and there are no systematic errors (Beck and Arnold, 
1977). In such situations, the smoothening of 'noise’ in the 
raw data can be accomplished by adopting a 'least square 
polynomial’ rather than an exact polynomial, provided the 
adopted polynomial closely approximates the true functional 
relation between the groundwater head and space coordinates in 
the given domain of space. The relevance of smoothening 
processes can be directly Inferred from the smooth water table 
or piezometric surfaces encountered in physical situations of 
alluvial aquifers ( Neuman, 1973; Sagar et al , 1975). 

Kriging, a spatial interpolation scheme, was used to 
obtain the head variations in the flow domain (Yeh, Yoon and 
Lee, 1983) as a presampling filter to reconstruct the head 
distributions over the spatial computational grids (or nodes). 
Kriging, when used with direct approach (Sagar, 1975) for 
parameter identification, replaces repetitive numerical 
simulations of the governing equations in the optimization 
approach . 

Kitanldis and Vomvorls (1983) applied Kriging to provide 
minimum variance and unbiased point estimates. The statistical 
approach developed by them has the significant advantage in 
that it does not attempt to restore details which can not be 
actually extracted from available data. 

3.4 PROCEDURE ADOPTED 

The locations for the available groundwater data 
generally would not coincide with all the grid points . So it 



is necessary to interpolate the missing groundwater heads at 
the nodal points of the computational grid with the help of 
least square approximation followed by statistical tests. The 
coordinates (x,y) of the observation points and head values at 
those points constitute the input data. Before going for 
interpolation it is mandatory to normalize observation points. 


3.4.1 Least Square Polynomial Approximation 

The spatial and temporal variation of groundwater head in 
an unconflned aquifer is governed by an equation of the 
diffusion type. Therefore, one may consider the observed head 
values available at finite number of points in space and time 
to be made up of two components i.e. 
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~ ^ (P ^ E 1 = 1,2, ...n 

k= 1 , 2 , . . . m 


(3.1) 


where k = observed head at i observation point (P i ,q i ) 

and at k th time, 

h k = the true v& l u ® head at (Pj.q^) and at k th 

discrete time, 

(p,q) = function of space coordinates, 

p,q = solution of groundwater flow equation at k th 
discrete time, 

E - random error associated with observations 

and assumed to be statistically un-correlated, 
n = number of observation points, 

and a = number of time periods 


For real situations, it may not be possible to know the 
functions of h k (p,q) explicitly. However, it is known a priori 



that these functions are continuous in space. The continuity 
of these functions can be directly inferred from the 
continuity of water table or piezometric surfaces which do not 
show any spatial discontinuity in aquifers free of 
hydrogeological structures. Thus, the function h k (p,q) can be 
approximated by polynomials of spatial coordinates p and q, 
and by H k (p,q) (least square polynomial solution). By choosing 
an appropriate degree of the approximating polynomial, the 
difference between the true function and the approximating one 
can be restricted to a small but non-zero value i.e. 

|h k (p,q) - H k (p , q) | | < <5^ (3.2) 

k - 1 to m and <5 > 0, where <5 is the error 

e e 


For obtaining a least square solution, the total number 
of the coefficients in a chosen polynomial must be less than 
the total number of data points. The coefficients of the 
approximating polynomials are estimated employing the observed 
groundwater head data (fc^ k ), by the least square criterion 
which implies the estimation of the coefficients in the chosen 
equation such that the sum the squares of the deviations of 
the observed values from the numerically predicted ones is 
minimized i.e. 


minimize 



H. 
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The residues at the 


l t * 1 observation 


point at 


th 


time 


point are given by 
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•i = h i.k - ®k <vi } 

A complete estimation of errors (E) and one to one 
representation of (p,q) by (p,q) can be ensured only if 
the residues (e^ are Identical to the errors in the raw data 
(E). However, such a situation does not generally occur 
because of the following reasons : 

(1) The true form of the function h k (p,q) is not known 
and it is being approximated by H^(p,q). According to the 
Uleratraas theorem, the minimum difference between the two 
functions within the domain of the observation, will always be 
grater than zero. 

(2) Even if the true form of function h^Cp.q) were 
unknown, (e^) will be identical to (E) only if the latter are 
normally distributed. Because of the finite sample size, the 
frequency distribution of the errors may be approximately 
normally distributed. 

3.4.2 Degree Of Polynomial 

As the degree of polynomial is increased, it may get 
closer and closer to the true functional relation between the 
groundwater head and space coordinates in a given domain. 
However, beyond a certain limit, any further increase in the 
degree of polynomial may only Induce the approximating least 
square polynomial to conform to the 'noise’ pattern of the raw 
data Instead of providing a better approximation of the true 
groundwater head. A partial polynomial may provide a better 
approximation than a full polynomial. It is desirable to 



arrive at an optimal form of the polynomial i.e. a partial 
polynomial of minimum possible degree and terms. Thus, apart 
from providing a better approximation, it also minimizes the 
computational efforts. The optimal form of the approximating 
function can be arrived at by employing the statistical tests 
of significance (Daniel and Uood, 1971). The standard tests of 
significance used in the context of multiple regression, can 
be directly applicable to least square polynomials provided 
the 'independent variables’ (i.e. the terms containing the 
different exponents of the spatial coordinates) are not 
correlated to each other (Dooge, 1973). This is applicable if 
the spatial coordinates and hence the independent variables 
are known exactly (without any error) and are thus, 
deterministic in nature. 

The degree of polynomial can be decided on the basis of 
criterion on minimum standard error a. 

2 2 

s = £ e i / ( n " n P ) where n>np (3.5) 

i 1 

where 

e^ = residue at i** 1 observation point (Eqn 3.4), 

2 

£ e^ = sum of the squares of the residues (SSR), 
n = total number of the data points, 

(n - np) = the degree of freedom (DF) of the 
proposed least square polynomial, 
and np - number of coefficients in the polynomial. 

As the degree of the polynomial is Increased, there will 
be a reduction in SSR but at the same time DF will also 



bu 

decrease since the third degree and a fourth degree full 

polynonlals would have ten and fifteen terns respectively. The 

/ 

polynomial which gives minimum standard error (s) la 
tentatively chosen as the approximating function (Ralston, 
1965) . 

3.4.3 Truncation Of Polynomial 

For any given level of confidence, the confidence limit 
for each of the estimated coefficients can be defined. This 
criterion can be used to delineate the terms which contribute 
significantly towards the explanation of the power of 
polynomial model. 

The confidence limits of estimate for coefficient of the 
1*** term are given by ; 

bj ± t ( n - np , a/2 ) s.e (b^) (3-4) 

where t( n-np, a/2) is (1-a) percentile of t-distrlbution 
with (n-np) degree of freedom. 

The standard error of b^, s.e. (b^) ia given by 

s.e (t^) = a V c ii (3.7) 

where c^ is the 1 th diagonal element of the corrected 
sum of squares and product matrix. 

The statistics for testing (b^=0) is 


t = 


s.e (b A ) 


(3.8) 
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Null hypothesis is accepted If the computed value of 't’ 
is lower than the critical value of 't’ for the degree of 
freedom (n-np) at some chosen level of confidence (generally 
95 percent). Acceptance of null hypothesis implies that the 
i** 1 term does not have significant explaining power and may be 
dropped. Thus, a few terms of polynomial can be dropped by 
adopting the students 't’ test. The test is described by 

* - 

t = (3.9) 

a / V~n 

where x = mean of the sample, 

= hypothetical mean of the population 

(18 percent), 

n = number of observations, 
and s = standard deviations of the observations. 

By this test, a reduced model can be prepared. However, 
it has to be tested if the reduced model (RM) gives as good of 
fit as full model (FM). Coefficients of the polynomial with 
reduced number of terms and SSR are reestlmated. If is the 
number of terms dropped then, 'F’ statistic is estimated by 


[SSR (RM> - SSR (FM> ] /K d 
SSR CFM)/ Cn - np~ > 


(3.10) 


If the calculated F is large in comparison to thetabulated 
values of F with K . and (n-np) degrees of freedom at a 

U 

percentile 1 evel of confidence , the result is significant at 

I 

level ot i.e. the RH is unsatisfactory. It implies then, even 

though each of the terms dropped indlvidualiyt nWtBRAJRrt 
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significant explaining power, collectively these terns explain 
a significant part of the variation in the groundwater head. 

3.4.4 Acceptability of the Approxinat ion 

Subsequent to the estimation of coefficients, it is 
desirable to make an assessment of the goodness of fit. The 
most widely used index is the multiple correlation coefficient 
R, defined by 

2 SSR 

R - 1 - (3.11) 

h I.k - V 

* 

n h i k 

where h = £ — (3.12) 

i = l n 

2 

R represents the fraction of the initial variance 
explained by the model. The adequacy of the proposed fit to 
approximate the true function may be checked by examining the 
residues. The standard residue at 1 th station is given by 

es jL = e x / s (3.13) 

The standard residues have zero mean and unit standard 
deviation. In general when the model is correct, the standard 
residues fall between +2 and -2 and their plots against 
Independent variables l.e. x and y do not exhibit any trend. 
Even after satisfying this criterion for the standard 
residues, some anomalous values still exist at some space 
points due to the presence of hydraulic discontinuities in the 
aquifer. They originate from geological structures like 



faults, folds and dykes. As a result of such structures, the 
water table will show narked differences In its elevations. 
Such groundwater heads can not be adequately approximated by 
least square polynomials. The discontinuities as evidenced 
from the contour map can be related to causative features 
identified from geological maps or remote sensed data. In 
these situations, estimation of derivatives can be monitored 
on either side of the structure by creating a temporary 
boundary between two space points in the problem domain. 

If the least square polynomial provides a good fit to 
the observed data, the hydraulic continuity of the aquifer is 
established. The 'goodness' of least square fit is 
characterised by high correlation coefficient (low standard 
errors), random distribution of standard residues in space and 
absence of outliers. The orientation and distribution of 
outliers can give some idea about the location and delineation 
of discontinuities, if present. However, the occurrence of 
outliers without any systematic orientation in space can be 
due to some localized phenomena like Incomplete recuperation 
at the time of observation and the inadequate hydraulic 
connection between the well and the aquifer. 

3.S PROCESSSING OF THE WATER TABLE DATA IN THE STUDY AREA 

The water table data used in this study are collected 
from Andhra Pradesh State Groundwater Department, Hyderabad. 
This data pertaining to a period of ten years (June78 to 
Hay88) was at fifteen day Intervals. There are 69 observation 
wells In the study area (Fig. 3.1 and Table 3.1). The data is 

















Table 3.1 : Details of Observation Uells Chosen for Data 


in the Present Study 


SI .No. 

location 

depth below ground 
level (m) 

f ormatlon(a) 
tapping 

1. 

Yekalaskhana peta 

8.05 

granite 

2. 

Thalleballl 

7.35 

granite 

3. 

Ratnapuram 

11.00 

granite 

4. 

Vasanthapuram 

8.30 

granite 

5. 

Trlpuraram 

8.00 

granite 

6 

Anataglri 

9.60 

granite 

7. 

Barakathgudem 

7.75 

granite 

8 . 

Akupamula 

9.40 

granite 

9. 

Komar abanda 

8.20 

granite 

10 . 

Kahanpur 

8.20 

granite 

11 . 

Shant inagar 

7.15 

granite 

12 . 

Dacharam 

7.60 

granite 

13. 

Lakkavaram 

9.20 

granite 

14 . 

Bethavolu 

15.00 

granite 

15. 

Mukundapuram 

8.10 

granite 

16 . 

Kodad 

7.30 

granite 

17 . 

Tammaband 1 apal em 

7.80 

granite 

18. 

Dorakuntla 

6.80 

granite 

19. 

Somavaram 

8.75 

granite 

20. 

Somavaram Tandla 

7.65 

granite 

21 . 

Fathepur 

6.60 

granite 

22 . 

Ponugodu 

10.90 

granite 

23. 

Thallamalkapuram 

4.65 

granite 

24 . 

Chilkooru 

6.25 

granite 

25. 

Kothagudlbanda 

10.00 

granite 

26 . 

Kapugallu 

8.50 

granite 

27. 

Chi llepal 11 

6.52 

granite 

28. 

Neereddlcherla 

6.45 

granite 

29 . 

Gardepalll 

9.20 

granite 

30 . 

Keethavarlgudem 

10.00 

granite 

31 . 

Rayanagudem 

10.40 

granite 

32. 

Gopalapuram 

9.60 

granite 

33 . 

Boorugugadda 

8.60 

granite 

34 . 

Huzurnagar 

9.80 

granite 

35 . 

Togarai 

6.70 

granite 

36 . 

Dirsencherla 

8.30 

granite 

37 . 

Garlkuntla Tanda 

7.20 

granite 

38. 

Yelapalasingaram 

9.65 

granite 

39. 

Llngaglrl 

10.35 

granite 

40. 

Yerlyaram 

6.35 

granite 


. . . Contd . . . . 




Table. 3.1 : Contd... 


SI. No. 

location 

depth below 
level (n) 

ground fornatlon(s) 

tapping 

41 . 

Ganapavaram 

9.50 

granite 

42. 

Yellapuran 

8.60 

quart 2 arenlte 

43. 

Palakeedu 

8.25 

granite 

44 . 

Gudlguntlapalem 

8.60 

granite 

45. 

Anar ar an 

10.75 

granite 

46 . 

Srlnlvasapuran 

8.20 

granite 

47 . 

Anandnagar 

8.61 

granite 

48. 

Nedlacheruvu 

11.90 

granite 

49. 

Kandibanda 

9.85 

quartz arenlte 

50. 

Revoor 

6.55 

linestone 

51 . 

Dondapadu 

10.90 

limestone 

52. 

Konatikuntla 

10.45 

quartz arenlte 

53. 

Bothapalen 

7.80 

granite 

54 . 

Kanalacheruvu 

6.00 

granite 

55. 

A1 angapur an 

6.60 

granite 

56. 

Yethavakella 

8.50 

granite 

57. 

Alllpuran 

6.95 

granite 

58. 

Nattanpalli 

9.00 

quartz arenlte 

59. 

Choutapalll 

13.65 

granite 

60. 

Gudlnalkapuran 

16.90 

linestone 

61 . 

Zonapadu 

14.70 

limestone 

62. 

Gundlapadu 

7.20 

quartz arenlte 

63. 

Raghunadhapal en 

9.50 

linestone 

64 . 

Hal lar eddlgudem 

22.15 

linestone with shale 

65. 

Peddavedu 

11.00 

quartz arenlte 

66 . 

Shobhanadr iguden 

14.10 

limestone with shale 

67. 

Chltryala 

20.20 

limestone 

68. 

Nakkaguden 

13.85 

linestone with shale 

69. 

Tannvaran 

12.95 

linestone with shale 


Ufel 1 locations are indicated in Fig. 3.1 with the sane serial 


nunbers . 



divided into two slabs of five years each. The data from both 
the five year seta (1978 to 1983 and 1983 to 1988) were uaed 
to estimate aquifer parameters by inverse problem. These 
estimated parameter values in turn have been utilized in 
calculating the water table elevations for the same periods. 
The estimated water table elevations were compared with the 
observed values to ascertain the accuracy of the model. The 
procedure involved in water table data processing is described 
in the succeeding pages. 

Least square polynomials were generated to approximate 
the spatial variation of water table elevation in sixty 
sequential months from June 1978 to May 1983. By this 
approximation a coefficient matrix was calculated. This 
coefficient matrix was employed to estimate spatial and 
temporal derivatives of hydraulic head. The axes (I,J) of the 
coordinate system along with the origin are indicated in 
Fig. 3.1. 

Procedure for the estimation of the minimum standard 
error to determine the degree of polynomial has been described 
(Section 3.4.2). On the basis of this error as a first 
attempt , coefficients of the third degree were estimated. A 
full polynomial of third degree would have ten terms for each 
period, as expressed below. 

H k (p,q) = bjp 3 + b 2 q 3 + b 3 p 2 q + b 4 p q 2 + b & p 2 + b fi q 2 + 

b ? p q + b 8 p + b 9 q + b 1Q 


(3.14) 


Using the above polynomial, the summary statistics were 
calculated for all the sixty months. The space points which 
displayed a standard residue of ±2 are treated as outliers. 
The 't' values were computed using the equation 3.8. These 
computed values of ‘t* are compared with the critical values 
of 't* at 95 percent confidence level and at degree of the 
freedom equalling the number of data points retained minus the 
number of terms. It was found that values of computed 't* of 
almost all the terms other than the constant term were lower 
• them - than the corresponding critical 't* values. Depending 
upon this criteria, some of the space points were deleted and 
the coefficients were once aglan calculated. This procedure 
helped in increasing the value of the variance but the 
standard error is too large to be acceptable. Moreover, large 
standard error indicates presence of large residues. As a next 
step, the polynomial terms showing large residues are deleted 
together with the apace points and the coefficients are once 
again estimated. The output indicated the standard error to be 
within the range of ±2 without any outliers. Even the space 
points showed improvement with this partial polynomial. But to 
achieve the satisfactory values of above criterion, many space 
points were deleted and polynomial equation had to be 
truncated. 

To Include all the values of space points and to have a 
better approximation, it was necessary to go for a higher 
degree of polynomial. A fourth degree polynomial having 
fifteen terms was considered. The full polynomial of fourth 
degree is of the form 



O y 


H^Cp.q) = b^p*+ b2^^+ b^p^q + b^p q^ + b^p^q^+ b^p^+ b^q' 


b e p 2 q + b,p , 2 + b 10 p 2 t b^q 2 . b^p q ♦ 


b 13 p + ‘“k"’ * b 15 


(3.15) 


The coefficients of polynomials are calculated, using the 
same criterion as in the third degree case. This scheme is 
also tested for decreasing the number of outliers and 
minimizing the standard error. In this case too, some space 
points are deleted and some terms of polynomial are not 
considered after correlating their ' t ’ values with the 
critical *t* values at 95 percent confidence level. This 
truncated polynomial exhibited better spatial distribution and 
had lesser number of outliers than full polynomial. Uhen it is 
tested for F test, the results are encouraging and revealed 
that the dropping of some space points and/or the polynomial 
terms leads to an Improved approximation of spatial variation 
of water table data. 

As an example, the data for the period June 1978 to Hay 
1979 can be examined. These data at 69 observation wells were 
tested for the above procedure. The least square polynomials 
for the twelve months were calculated in the beginning with a 
third degree full polynomial with all the ten terms and at all 
the 69 space points. These coefficients are tested for summary 
statistics and the results are shown in Tables 3.2 and 3.3. 
The computed 't’ statistics for the month of February’79 are 
shown below as a typical example. 
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Table 3. 

2 : Deletion 

Scheme of 

Space 

Points for 

a Third degree 

Polynomial 

SI. No. 

month number of 

variance 

goodness 

of fit 


points 

deleted 

initial final 

initial 

final 

1 . 

Jun 78 

16 

0.9102 

0.9742 

81.55 

95.33 

2. 

Jul 78 

16 

0.9256 

0.9820 

81.32 

95.28 

3. 

Aug 78 

17 

0.9009 

0.9725 

81.54 

95.46 

4. 

Sep 78 

16 

0.8690 

0.9537 

82.24 

95.48 

5. 

Oct 78 

17 

0.9033 

0.9621 

81.18 

95.35 

6. 

Nov 78 

17 

0.9033 

0.9680 

82.84 

95.42 

7. 

Dec 78 

15 

0.8704 

0.9593 

82.14 

95.42 

8. 

Jan 79 

15 

0.8743 

0.9602 

82.27 

95.83 

9. 

Feb 79 

15 

0.8910 

0.9714 

82.04 

95.68 

10. 

Mar 79 

15 

0.8975 

0.9657 

82.05 

95.55 

11 . 

Apr 79 

15 

0.8975 

0.9648 

82.08 

95.57 

12. 

May 79 

15 

0.9067 

0.9733 

81.93 

95.23 


The common space points deleted are of serial numbers 
4,8,6,11,13,15,18,25, 26,27,35,38,39,55 and 67 on Fig. 3.1. 
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Table 3.3 : Truncation of Polynomial Terms for a Third Degree 

Polynomial 


SI .No. 

month 

number of 

terms 

del et ed 

variance 

standard error (m) 

initial final 

initial 

final 

1 . 

Jun 

78 

3 

0.9102 

0.9742 

0.9032 

0.7791 

2. 

Jul 

78 

3 

0.9256 

0.9820 

0.9222 

0.7216 

3. 

Aug 

78 

2 

0.9009 

0.9725 

0.8937 

0.6328 

4. 

Sep 

78 

2 

0.8690 

0.9537 

0.9072 

0.6219 

5. 

Oct 

78 

3 

0.9033 

0.9621 

0.8873 

0.6855 

6 . 

Nov 

78 

2 

0.9033 

0.9680 

0.9050 

0.7234 

7. 

Dec 

78 

1 

r 

0.8704 

0.9593 

0.8593 

0.5866 

8 . 

Jan 

79 

2 

0.8743 

0.9602 

0.9769 

0.6838 

9. 

Feb 

79 

3 

0.8910 

0.9714 

0.9012 

0.7243 

10. 

Mar 

79 

3 

0.8975 

0.9657 

0.9312 

0.7465 

11. 

Apr 

79 

2 

0.8975 

0.9648 

0.9060 

0.6791 

12. 

May 

79 

2 

0.9067 

0.9733 

0.8955 

0.6954 


2 2 

The common terms deleted are p q, pq and pq. 


Term 

t value 

3 


P 

-0.94 

3 

q 

-1.20 

2 

p q 

0.25 

2 

pq 

-0.34 

2 

p 

1.25 

2 

q 

-1.17 

pq 

-0.45 

p 

0.90 

q 

-2.36 

constant 

30.50 


The 

standard error Is 0.9012. 






The variance value Is about 0 

.8910 

and the 

goodness 

of 

fit 

Is estimated as 82.038. As the 

space 

points 

4/8, 6, 

11, 

13, 

15, 18, 25, 26, 27, 35, 39, 55 

and 

67 have displayed 

deviation and presence of outliers. 

they 

were deleted and 

the 

coefflcents were computed again. The statistics 

of fit for 

the 

rest 

of the points Is as follows: 






Term 

t 

value 




3 

P 


1.22 




3 

q 


2.56 




2 

p q 


0.48 




2 

pq 


-0.34 




2 

P 


1.53 




2 

q 


-2.23 




pq 


0.89 




p 


-1.54 




q 


-3.45 




constant 


564.23 





A.a per the criterion, the ploynomlal t erne showing 

2 2 

residuals of ±1 i.e. p q, pq and pq are deleted. After this 

truncation, the coefficienta and the summary etatiatica were 

once again calculated. Uhile the standard error of partial 

polynomial fit worked out to be of a low value (0.7243 m) the 

variance value rose to 0.9714. The goodness of fit has also 

increased to 95.677. But to attain these results, 15 space 

points and three terms had to be curtailed. It may be noticed 

that in each month’s data, large number of reslduls (about 15 

points) are present (Table 3.2). Thus the inadequacey of third 

degree polynomial in approximating the spatial variation of 

water table during these months is established. This 

neccesslated the recomputation of the least square polynomials 

for all the months using the fourth degree polynomial. This 

polynomial was truncated on the lines as described earlier. 

The space points truncated are 4 , 13, 35 and 67. The final 

fourth degree partial polynomial has got 13 terms and 65 space 

3 2 

ponlts, after deleting the terms p q and p depending on the 
criterion adopted for the removal of residuals. 

This partial fourth degree poynomial equation was used to 
approximate the spatial variation of water table head values 
at all 69 points in the study area for the 120 months period. 
The head values are then interpolated using the equation 

h k (p i ’ q i> ~ H k (p ’ q) (3.16) 

The head values thus processed are stored for later use 
to assist in the calculation of the second spatial derivatives 
and the first temporal derivatives. Uith help of these 



derivatives, the inflows and outflows across the boundary and 
the changes in storage as well as the aquifer parameters are 
estimated, the processes for which are described in Chapters 4 
and 5 . 

The computer program to estimate the least square 
polynomials and to calculate groundwater system elements 
(Chapter 3 ) has been developed by the author in the course of 
present study. The algorithm for the same is presented in 
Appendix A. 

3. 6 CONTOURING AND TREND SURFACE PLOTTING 

Contour Haps 

Construction of contour maps of spatially distributed 
data presents problems similar to those in the creation of 
equally spaced data from the observations located at irregular 
intervals along a line, since the first step In contouring 
usually is to produce a regular mesh or grid of control points 
(Davis , 1973) . The regular grid may be produced in a variety of 
ways, ranging form estimates at grid points based on nearest 
observation to estimates derived from trend surfaces of all 
observations. 

A set of values to be contoured is entered into the 
computer as (3 X n) array, where each data entry corresponds 
to X. coordinate, X. coordinate, and Y or dependant variable 
which will be contoured. Thus, the Inputs to the contour 
program are given in terms of an east-west coordinate (X^), a 
north-south coordinate (X 2 ) and elevation of water table above 
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datum plane (Y) for each of the observation points which are 
expressed by numerical coordinate system. Then, a rectangular 
grid of points is created which will control the computing 
process. This grid will become a rectangular mesh of equally 
spaced points estimated from the data. The grid is determined 
by the equation 


n 


E a,/D ik ) 


Y k = 


i = l 


n 


E Cl/D ) 
i = l 


C3-17) 


* A.C 

where Y is estimated elevation at k grid point, n is 
the number of observation points and is the distance from 
observation point i to grid point k. 


D 


ik 


V cx 


lk 


Y -’> + (Xzk'^i 5 


li 


. (3.18) 


where X lk and are grid points corresponding to X^ 
and X 2i observation points. 

In ordinary contouring routine, after establishing grid 
points and corresponding values the next step would be to 
interpolate between adjacent grid points and determine the 
coordinates of any contour line that may pass between two 
points and this is done by linear interpolation. A simple 
algorithm for printing contour maps on a line printer can be 
adopted or with the help of a graphical package, contours can 
be drawn using a dot-matrix printer. 
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The 

contour map 

is 

helpful 

in 

finding 

out 

water 

table 

position 

at each grid 

point . 

Pronounced 

changes 

in the 

hydraulic 

continuity 

as 

evidenced 

on such 

a 

map 

can be 


attributed to causative geological features, if present, or to 
the erroneous interpolation during least square polynomial 
approximation. 

Trend Surface Plotting 

A trend surface is a surface generated by a linear 
function of the geographical coordinates of a set of 
observations, so constructed that the squared deviations from 
the trend are minimized (Davis, 1973). The trend is a linear 
function. That 1s, it has a form Y = bX^ + bX 2 + ...., where 
the b’s are coefficients and the X's are some combination of 
the geographic coordinates. The equation would yield values of 

A 

the Y which are trend components of an observation such as the 
groundwater head. The specific linear function chosen for 
trend must minimise the squared deviations from the trend. A 
linear trend surface is represented by an equation of the type 

* * V b lV b 2 X 2 <31,) 

A geologic observation, Y, may be regarded as a linear 

function of some constant value (b ) which is defined as 

« 

/ 

b = f (b. ,b.) + n (3.20) 

O JL 4w 

where b., and b_ are east-vest and north-south components 
1 z 

and n is the number of observations. Since three unknowns are 
involved in this equation, three simultaneous equations are 
needed to solve the same. 
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E Y = t£+ x 1+ b 2 E x 2 

E XjY = b o E Xj+ *> ± E X* + b 2 E X jL X 2 

E x 2 y = b o E x 2 + bjE b 2 E x 2 

The above equations can be rewritten 


(3.21) 


in the matrix form as 


”n 

E x ! E X 2 1 


b ol 


■ ey 

E 

E x i E x i x 2 


b l 

= 

EX X Y 

X X 2 

E X 1 X 2 EX 2 2 J 


b 2 J 


. X X 2 Y . 


A least square line may be expanded to a second degree 
curve (parabola) by adding a squared term to the linear 
equation. 

Y = b 0 * bjX 1+ b 2 *\ (3.23) 

An equivalent expression of (Eqn.3.23) gives a 
second-degree trend surface 

Y = b Q + b 1 X 1 + b 2 X 2 + b 3 xj+ b 4 X* + 

(3.24) 

This procedure can be extended for higher degree of trend 
surfaces, with more unknown coefficients. 


A trend surface program consists of three parts: 

(1) A routine to generate the matrix of sums and cross 
products . 

(2) A simultaneous equation solver or a matrix inverter, 
and (3) A plot program. 

Grid points (X^, X 2 values) .dimensions of the map and 
scaling factor (D^) can be calculated as In contouring 



program. But the observation values (Y) are employed, after 
fitting them to a second or higher degree of trend surface 
equation, to generate a trend surface map. The goodness of fit 
of a trend surface may be tested statistically, by comparing 
the variance due to regression or trend to the variance due to 
deviations from the trend. The/sontour plots are prepared with 
the help of a graphical package (PLOT 88) in MSFORTRAN and the 
trend surface diagrams are generated on IBM PC/AT with the 
help of a EGA driver on the basis of the program written in C 
language as developed by the author. 

3.7 CONFIGURATION OF THE WATER TABLE IN THE STUDY AREA 

Uater table contour maps for the months of April, August 
and December of 1981 are presented (Figs. 3.2, 3.3 and 3.4.) to 
indicate typical pattern for one year in the study area. Uhile 
there is a variation of the head values in the three cases at 
different points, the patterns are in general similar in all 
the three cases. Several distinct trends can be noticed. A 
general gradient of the water table head from north to south 
occurs with variation from over 150m (above MSL) to 56m (above 
NSL) in the southeastern corner of the area. 

Effect of geology on the water table configuration is 
evidenced on overlapping the water table map on a geological 
map. For example in the area where dykes are present, a 
modification of the slope is pronounced as in the areas south 
of Ponugodu, east of Kamalacheruvu, north west of Llngagirl 
and south west of Ponugodu (Fig. 3. 2). A marked change in the 
head values in the vicinity of the dykes is noticeable in the 



granitic terrain within these areas. J 

Effect of lineaments on the water table head 
configuration has also been evidenced. In the granitic terrain 
west of Bethavolu, marked change in the water table contour 
plots can be seen in the region where the lineament la 
located (Fig. 3. 2). Similarly within the limestone terrain a 
NE-SU lineament in the area north of Chltryala appearea to 
have influnced the contours in widening their spacing 
indicative of increased permeability. 

Typical trend surface plots for water table head 
variation at 4 years intervals are indicated in Figs. 3. 5, 3.6 
and 3.7. It is observed that in all the three cases, the 
trends are similar. As seen in the contour maps, although a 
change in the head values exist at different periods, while 
the present vertical scale of the trend surface plots the 
differences are not pronouncedly brought out in the figures. To 
provide an idea of topography in the area, the mesh for 
topography has also been incorporated (white grids in 
Fig. 3. 5). 

Thus, in the present chapter reconstruction of water 
table head data has been carried out by minimizing sum of the 
squares of the residual errors. A functional approximation was 
considered with necessary statistical tests to estimate the 


least square 

polynomials . 

The 

data 

were processed 

and 

the 

coefficients 

were calculated. 

These 

data were 

later 

used 

to 

estimate the 

groundwater 

system 

el ements 

(Chapter 

4 ) 

and 


aquifer parameters (Chapter 5). 
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3.5 


Uater Table Head Configuration December, 1979 


/ -> 



Water Table Head Configuration December, 1983 
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3.7 


Uater Table Head 


Configuration 


December , 


1987 


CHAPTER 4 


ESTIMATION OF GROUNDWATER RECHARGE AND DISCHARGE COMPONENTS 

4. 1 INTRODUCTION 

For proper evaluation of aquifer response to any given 
excitation, the various recharge and discharge elements are to 
be estimated. These elements are from the surface as well as 
groundwater source. Rainfall, canal seepage, evaporation and 
irrigation effect are from surface source while the 
groundwater system elements include flow across the boundary, 
storage and well pumpage. The algebraic sum of these values 
form an Important input in addition to the water table 
elevation data in groundwater modelling. The groundwater head 
values obtained by processing (Chapter 3) are used for the 
estimation of flow across the boundaries and storage. In this 
chapter, procedures for the estimation of all the system 
elements are presented and the same are estimated for the 
study area. 

4.2 ELEMENTS OF GROUNDWATER BALANCE 

The hydrologic equation based on the law of conservation 
of matter, as applied to hydrologic cycle defines the total 



water balance. Groundwater balance deala with aspects of 
balancing various components of groundwater supply (recharge) 
and disposal (discharge) with storage changes in the ground 
water reservoir. The items of recharge and discharge are as 
follows : 

Components of recharge 

1. Precipitation, infiltrating to the aquifer 

2. Natural recharge from streams, lakes and ponds 

3. Groundwater inflow 

4. Artificial recharge from canals and irrigation 
Components of discharge 

1. Evaporation from the shallow aquifers and 
transpiration from plants. 

2. Natural discharge by effluent seepage, spring flow 
and base flow to streams and lakes. 

3. Groundwater outflow . 

4. Artificial discharge by pumping or through drains. 

The groundwater balance of a basin or an area, for an 

Inventory period may be expressed as 

Groundwater Inflow - groundwater outflow = 

change in groundwater storage. 

4.2.1 Estimation of Recharge Components 

The methods adopted for estimation of recharge take into 
account intensity and duration of rainfall, evapotrans- 
plratlon, soil moisture, runoff, infiltration capacities of 



soils, storage characteristics of aquifers and water level 
£ 1 tic tuat Iona and movement of ground water. The various 
recharge components to be estimated and the methods employed 
are discussed in the following paragraphs. 

Infiltration occurring at the land surface can be 
estimated by the soil moisture balance approach. The soil- 
moisture balance for any time interval can be expressed as 


P s AE + I + R +AS 

in 

where 

P = rainfall, 

AE = actual evapotranspiration , 

AS ra = change in soil moisture storage, 
I = infiltration, 
and R = surface runoff 


C4.1) 


Soil-moisture budgeting, taking into account evapo- 
transpi rational abstraction from precipitation, provides a 
measure of moisture available for runoff and infiltration. 
This can be done by Thornthwaite’ s formula. In this method 
measurements of field capacity and wilting point are to be 
measured from the available moisture down to the root zone. 

For long duration of rainfall and under conditions of 
excess of Infiltration capacity, the hydrologic equation can 
be expressed in the following form provided, the losses due to 
evapotranspiration are taken negleglble. 
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P = R + Up (4.2) 

where 

P = rainfall, 

R = surface runoff, 

and Up= recharge by infiltration from rainfall. 


The rainfall infiltration factor can be determined for 
catchments with varying terrain, soil and other 
characteristics . The following equation can be used to 
determine the rainfall-infiltration factor 


" <*bf + 


GU + 
e 


V 


+ U - U 
as ag 


- S ) 
P 


(4.3) 


where 


I = rainfall infiltration factor, 

P 

R, = base flow, 

DI 

GU = groundwater extraction, 

E t = evapotransplratlon from groundwater, 

E - evaporation from base flow, 

U = recharge due to infiltration of applied surface 
as 

water , 

U = recharge by return circulation from applied 
ag 

groundwater, 

and S = stream-bed percolation, 
p 
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Uater-level hydrographs can be decomposed to obtain 
fluctuations in response to rainfall and other inputs of water 
to the aquifer. The rise of water level in an aquifer 
represents the net response to a process of simultaneous 
drainage or discharge from and recharge to the aquifer. In 
order to find the actual recharge, the water-level recession 
has to be extended till the period of maximum rise of water 
level. If the recharge is discontinuous with large time gaps, 
the water level will decline between two recharge periods. The 
total recharge will be the sum of individual rises. The water 
table fluctuation observed in response to rainfall in a single 
year may not be representative of average long term 
fluctuations and recharge. The average rate of groundwater 
recharge to water table aquifer can be estimated from the 
shape of water table of an aquifer bounded on two sides by 
parallel streams (canals) of infinite length. If the aquifer 
is homogeneous and isotropic and is recharged at a constant 
rate of accretion in space and time, the rate of recharge (U) 
can be written in the form 

U = - 2 (4.4) 

a.T3*io-»(“ -Ik ) 

where 

d = rate of recharge (cm/year), 

a = distance (m) from stream to the groundwater divide, 
x = distance (m) from the stream to an observation well. 



hp= elevation (m) of the water table at the observation well 

with respect to mean stream level, 

2 

and T - transmissivity (m /day). 

4.2.2 Recharge from Canals and Streams 

Recharge through percolation from canals and streams 
depends on the infiltration capacity of the canal and stream 
bed and sides as well as subsurface lithology, extent of 
wetted perimeter, length of canal or stream, discharge, 
sediment load, physical and chemical properties of water, and 
relative position of the bed with respect to the water table. 
As in the case of streams, canals may be Influent or effluent 
and it is essential to establish the relation before and 
during the studies. Recharge rates may decline over the years 
due to water-logging or clogging of pores of the bed material. 

For quantitative estimation of percolation rates, a 
stretch of canal or stream with no diversion of water is 
generally selected. In case of canals, the seepage losses are 
expressed usually in m 3 /s per million sq m of wetted surface 

as 

Oj, ~ °2 (4-5) 

A 

where 

I = seepage losses in cumecs per million sq.m, 

Q = discharge at upstream point, cumecs, 

1 

0 = discharge at downstream point, cumecs, 

2 

and A = area of wetted surface (sq.km), (length of canal X 



average of perimeters between two discharge points). 

Recharge from canals and streams that are in direct 
hydraulic connection with a phreatic aquifer underlain by a 
horizontal impermeable layer at shallow depth can be 
determined by Darcy’s equation, provided the flow satisfies 
Dupuit assumptions. 



Uhere h and h. are saturated thickness values of the 
s 1 

aquifer over a distance L at the stream or canal. For 
calculating the area of flow cross-section, the average of the 
saturated thickness (h + h. ) is usuaslly taken. 

S 1 

4.2.3 Estimation of Discharge Components 

Once the total groundwater recharge has been estimated, 
the groundwater losses have to be accounted for a quantitative 
estimation of the exploitable surplus. 

For the estimation of subsurface outflow of groundwater, 
contour maps of groundwater surface prepared on the basis of 
the water table level data of wells located both within and 
outside the section delimiting the basin outlet are used. 

The quantity of groundwater discharge appearing as base 
flow can be measured by stream gauging at the basin outlet. 
Evapotransplratlon from fully saturated soils with vegetation 
equals the potential evapotransplratlon rate. Evaporation loss 
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from soil surface forms one of the major components of 
groundwater losses and Is pronounced in the water-logged 
canal command areas. Draft from individual wells may vary 
widely depending on the yield, type of well, well design, 
depth of water level, method of lift, crops grown etc. 

1. 3 RECHARGE CALCULATION 


In general, the rainfall recharge and groundwater pumpage 
are the dominant factors influenceing the recharge estimation. 
However, in canal command areas seepage from canals 
contributes significantly to the net recharge along with other 
components of Q (total algebraic sum of inflows and outflows). 
Thus, for Q in 1 space and K time period can be written as 


r ik 


= R - U.. + C.. 

Ik lk ik 


+ X 


ik 


(4.7) 


where R ifc is the rainfall recharge, is the effective 
wlthdrawl due to pumpage (with adjustment to the return 
flows), is recharge from canals through seepage and X ifc is 
algebraic sum of all other X ik recharge components of 0 i j c » 
with the recharge taken as positive and discharge taken as 
negative values respectively. 

In tropical regions, the recharge from rainfall is an 
important component of groundwater balance. The available 
rainfall records are usually adequate to estimate the spatial 
and temporal distribution of rainfall. 



85 

In the present study, the following equation proposed by 
Kashyap (1981) was used to calculate recharge elements of the 
aquifers. 

R ik " f r [ P i,k’ P i,k-1’ P i ,k-ne’ a ) (4-8) 


where ^ is the rainfall at i** 1 space point during k th 
period, me is the number of periods preceding rainfall 
affecting the recharge in the current period and oi i is a row 
matrix containing number of parameters. In the procedure for 
soil moisture budgeting (Eqn.2.7), the parameter determined 
will be the field capacity. In the absence of the availability 
of the bulk data, the simple functional form that can be used 
for the estimation of the built-in parameters is 


R.. = K P.. 

Ik r lk Ik 


(4-9) 


where K is the recharge coefficient and P.. is the 
r ik iK 

precipitation. It varies from basin to basin, and for any 


given basin, a temporal variability of K 


exists. For 


ik 


shorter time periods, the antecedent soil moisture conditions 
can be incorporated through a parameter threshold rainfall 
(«*) which is the minimum rainfall required in one time period 
to fill the soil moisture deficiency. Thus, rainfall recharge 


relation can be written as 
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R ik = °1 


= 0 


if p i,k-i< 


2 P i.k + (' - “ 2 ] P l.k -1 - “3 
lf “3 > P l,k * ( 1 ' “2 ) P i.k -1 » n<1 P i,k -1 < <4 


Nou *lk = “l 


“2 P i,k + ( 1_ “ 2 ] P I,k-l 


3 

(4.10) 
i 


i£ P i,k-1 * "3 


(4.11) 


where ot £ and a* are recharge parameters of the i** 1 
point. The a* corresponds to the recharge coefficient and a* 
accounts for the delayed aquifer response. The approach for 
incorporating the rainfall recharge in the Inverse problem 
formulation not only provides consistent estimates of 
historical rainfall recharge, but also effectively calibrates 
the locally acceptable functional form of the rainfall 
recharge relations (Kashyap, 1981). 

The recharge due to canal seepage can be calculated from 
the expression (modified form of Eqn.4.3) 


C lk - ! c + l r ♦ "l ♦ *w - D 0 


.- (4.12) 


where I c is recharge due to percolation from main canal 
and branches, I is the recharge due to percolation from 
reaches, I A is the percolation from canal irrigated areas, 1^ 
is percolation water released from areas Irrigated by wells 
and D q Is the discharge from canals, branches and reaches 
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across the boundary of basin to outside of system under study. 
Thus, total loss of groundwater to the rivers as subsurface 
outflow from the area is also included in this component Dq. 

<4. 4 ESTIMATION OF THE GROUNDWATER SYSTEM ELEMENTS 

As indicted earlier, (Chapter 3), the approximating least 
square function is employed for obtaining system elements 
relating to the groundwater studies. 

4.4.1 Groundwater Storage 

The groundwater storage in k time can be estimated by 
integrating the function H k (p,q) over the space defined by the 
study area. Thus, the groundwater storage SG above a datum 
plane ZD is given by 

SG fc = J J s[ H k (p , q) - ZDj dp dq (4.13) 

AG 

where AG is the study area. The integration over an 

irregularly shaped area can be carried out by dividing the 

/ 

area into finite number of strips extending along one of 

two directions p or q. The groundwater storage in the 

atrip (Fig. 4.1) is given by 

p? m+Aq 
i Q i 

H k (p.q) - ZD j dp dq 

p=p--q=q i 



the 

,th 


(4.14) 





(4.15) 


Total groundwater storage is given by 

ns 

SG = Z AS. 
i= 1 , n 1 

where ns is the number of strips. The function [ (p,q) 

ZD] can be integrated analytically. Thus, AS^ can be 
evaluated directly if the storativity £ is uniform in i th 

a. L 

strip. If S is varying in i strip, then SG may have to be 
evaluated by the numerical solution of equation (4.15). 


4.4.2 Boundary Recharge 

Total subsurface inflow rate (1^) across the boundary of 
the an area at k ttx time point is written as 


I 


k 



a 

Sn c 


H k (p,q) T dc 


(4.16) 


c 

where, c is boundary contour and n is its inner normal. 
The integration can be performed by dividing the boundary 
contour into finite number of elementary lengths (Fig. 4. 2). 
The subsurface flow over i** 1 elementary length is given by 

* = —4- H * <>vv AL i T i 


where p^,q^ are the two coordinates of the central point 
of I th elementary length and T^ is the average transmissivity 
at its elementary length (AL^). 



Orthogonal axis(q ) 


FIG. 4- 



Orthogonal axls(p) 


ESTIMATION OF LATERAL FLOWS ACROSS 
THE BOUNDARY 
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ne 

I k - E * I, (4.18) 

where ne la the total number of elementary lengths. In 

c 

equation 4.18, the hydraulic gradient in the direction of n 
can be evaluated directly, by knowing the average orientation 
of the elementary length with respect to p and q directions. 

In the present study, the polynomial coefficients are 
calculated first and and then groundwater storage and boundary 
recharge are estimated using the equations 4.14 to 4.18. 

4.5 ESTIMATION OF RECHARGE COMPONENTS IN THE STUDY AREA 


For the study area, the recharge has been estimated on 
the basis of data pertaining to rainfall and canal seepage, 
and discharge values from well draft and evaporation factors. 
The total algebraic sum of inflows and outflows (Q ik ) at * 
space point and k th period is expressed as 


'ik 


B — U + 

K ik ik 


ik 


X 


ik 


(4.19) 


The matrix R ik indicating recharge due to rainfall is 
calculated by the equation 4.19. In this equation, P lk is the 
precipitation during k th time period and at 1 th space point, 
and the data for the same were collected from Andhra Pradesh 
State Groundwater Department. There are only two rain gauge 



estimates. The matrix U ik (in Eqn.4.19) representing effective 

withdraw 1 due to pumpage. The draft from each individual well 

lS calculated using the norms prescribed by Andhra Pradesh 

Sate ground Uater Department. The discharge values considered 

from a tube well and from a dug well are 35000m 3 /day and 
3 

20000m /day respectively (including the effect of other 
pumping wells which are not observation wells). The effective 
withdrawl values due to pumpage (U.^) were calcq^flated using 
the above norms at 69 observation wells and for 120 months. 

Values of for the month of Harch , 81 are presented in 

Table 4.1. 

The term in Eqn.4.19 represents recharge through 

seepage from canals and reaches. Estimates of total seepage 
were obtained taking into consideration conveyance losses for 
the main canal and branch canals as 1.858 million cubic meters 
and 1.1148 thousand cubic meters respectively, based on the 
recommendation of George Smout , Uorld Bank expert (A.P.S.G. 
U.D. Report, 1981). The actual canal discharge values (in 

3 

m /month) for each canal have been calculated on the basis of 

the total discharge data (as obtained from Andhra Pradesh 

Irrigation Department, Canals Office of Nagarjuna Sagar Left 

t 

Canal Command area, Khammam, A.P) and the total number of 
running days. The recharge due to canal seepage for all the 
120 months at 69 observation points were estimated depending 
on their distance from the canals and/or reaches and and the 
polygonal area of Influence for the well. Estimates of C lk for 
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the month of March, '81 are also presented (Table 4.1). 

All the other recharge and discharge components are taken 
care of by the term (Eqn.4.19) where losses due to 
evaporation and irrigation effects on canal recharge etc., are 
included in the matrix. The value for the evaporation, taken 
as 0 . 987E-07m/sec , is based on the monthly evaporation data 
obtained from A.P.S.G.U. Department. The transpiration effect 
was not seperately considered. The irrigation effect is taken 
as 5.0* of discharge running in the canals and it is equally 
distributed to all the observation points. These values at all 
69 observation points and for all the 120 time periods are 
also calculated. The values for a typical month are shown in 
Table 4.1. 

The calculated recharge and discharge values from various 
factors such as rainfall, canal seepage, irrigation effect, 
groundwater wlthdrawl and evaporation are stored for use in 
the estimation of aquifer parameters (Chapter 5). 
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Table 4.1: Recharge Values for a Typical Month (March, 1981) 


(in thousand cubic meters) 


Observation 

R 

C 

U 

X 

0 

well No. 

1 

2 

3 

4 

l+2-(3+4 ) 

1 . 

55.58 

184.09 

20 

25 

194.67 

2. 

66.69 

184.09 

35 

25 

190.78 

3. 

116.71 

479.68 

20 

25 

151.39 

4. 

63.92 

479.68 

20 

25 

498.60 

5 . 

247.33 

783.06 

20 

25 

985.39 

6 . 

77.01 

303.37 

20 

25 

326.18 

7. 

61.14 

11774.42 

20 

25 

11790.56 

8. 

86.14 

184.09 

35 

25 

210.24 

9. 

63.92 

184.09 

20 

25 

203.01 

10. 

61.14 

479.68 

20 

25 

495.82 

11. 

18.53 

303.37 

20 

25 

276.90 

12. 

17.37 

274.85 

20 

25 

247.22 

13. 

69.47 

36.30 

35 

25 

45.77 

14. 

80.58 

36.30 

20 

25 

71.89 

15. 

111.16 

11738.17 

20 

25 

11804.33 

16. 

50.02 

184.09 

35 

25 

174.11 

17 . 

66.69 

479.68 

20 

25 

541.38 

18. 

34.09 

479.68 

20 

25 

458.77 

19. 

8.68 

1685.39 

20 

25 

1649.07 

20. 

11.58 

324.11 

20 

25 

290.69 

21 . 

20.26 

342.26 

20 

25 

317.52 

22 . 

24.90 

67.46 

20 

25 

47.36 

23. 

23.16 

95.53 

20 

25 

73.69 

24 . 

37.05 

11738.17 

20 

25 

11730.22 

25 . 

19.11 

11738.17 

35 

25 

11697.28 

26 . 

31.84 

11738.17 

20 

25 

11725.01 

27. 

4.05 

11738.17 

20 

25 

11697.22 

28. 

11.58 

1037.17 

20 

25 

1003.75 

29. 

22.00 

560.07 

20 

25 

537.07 

30. 

20.26 

95.53 

20 

25 

70.79 

31 . 

2.89 

95.53 

35 

25 

38.42 

32. 

5.21 

95.53 

20 

25 

55.74 

33. 

8.68 

11738.17 

20 

25 

11701.85 

34 . 

8.68 

11738.17 

20 

25 

11701.85 

35. 

21.42 

11738.17 

35 

25 

11699.59 

36 . 

21.42 

342.26 

20 

25 

318.68 

37 . 

32.42 

313.74 

20 

25 

301.16 

38. 

24.32 

344.86 

20 

25 

344.18 

39 . 

12.74 

11738.17 

35 

25 

11690.91 

40. 

9.84 

11738.17 

20 

25 

11703.01 


Contd 


Table 4.1 Contd. . . . 


96 


Observation 

R 

C 

u 

X 

o 

veil No. 

1 

2 

3 

4 

l+2-(3+4) 

41 . 

20.84 

11738.17 

35 

25 

11699.01 

42 . 

5.79 

342.26 

20 

25 

303.05 

43 . 

12.74 

342.26 ^ 

20 

25 

310.00 

44 . 

16.21 

342.26 

20 

25 

313.47 

45. 

19.68 

95.53 

20 

25 

115.21 

46 . 

12.16 

440.80 

20 

25 

407.96 

47 . 

17.95 

11738.17 

20 

25 

11711.12 

48. 

24.32 

11738.17 

35 

25 

11702.49 

49 . 

30.69 

11738.17 

20 

25 

11723.86 

50. 

42.26 

11738.17 

20 

25 

11735.43 

51 . 

37.05 

11738.17 

35 

25 

11715.22 

52 . 

13.89 

342.26 

20 

25 

311.16 

53. 

9.84 

342.26 

20 

25 

307.10 

54 . 

13.32 

342.26 

20 

25 

310.58 

55 . 

11.58 

342.26 

20 

25 

308.84 

56 . 

20.84 

11738.17 

20 

25 

11714.01 

57 . 

22.00 

11738.17 

20 

25 

11715.17 

58. 

15.63 

11738.17 

35 

25 

11693.80 

59 . 

21.42 

11738.17 

35 

25 

11699.59 

60 . 

30.11 

11738.17 

35 

25 

11708.28 

61 

39.95 

342.26 

35 

25 

322.21 

62 . 

23.74 

324.26 

20 

25 

302.99 

63. 

47.48 

11738.17 

20 

25 

11740.68 

64 . 

32.42 

11738.17 

35 

25 

11710.59 

65. 

63.69 

11738.17 

35 

25 

11742.86 

66 . 

41.69 

11738.17 

20 

25 

11738.85 

67 . 

23.74 

11738.17 

20 

25 

11716.91 

68. 

21.42 

11738.17 

20 

25 

11714.59 

69 . 

24.89 

11738.17 

35 

25 

11703.07 

R r Recharge 

due to rainfall 

C : 

Recharge 

due to canal 





seepage 


U : Discharge 

due to 

well draft 

X : 

Discharge 

due to evaporat 


and irrigation effect 


Q : Net inflows and outflows 


CHAPTER 5 


ESTIMATION OF AQUIFER PARAMETERS BY INVERSE PROBLEM 

B.l INTRODUCTION 

By known transmissivity (T) and storativity (S) , spat ial ly 
distributed over the problem domain, estimation of water 
levels by simulating the aquifer response to a deterministic 
pattern of ground water withdrawals and recharge is known as 
'direct problem* in groundwater hydrology. Pump test, a 
popular means of estimating aquifer parameters (T and S), 
Involves generating the aquifer response to the pumping in a 
single well. In the direct problem, T and S values, as 
obtained from pump tests are assigned to nodal points in a 
distributed model. However, these values are approximate as 
the pump test data analysis is based on several assumptions 
such as isotropic aquifer with a fully-penetrating well of 
infinitesimal size. Further spatial limitation of test pump 
test sites on a regional basis exisists due to the high cost 
of operations. 

On the other hand, the 'inverse problem* is an approach to 
estimate T and S by employing historical data of aquifer 
response and corresponding aquifer excitations. The aquifer 


/ 
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excitation can be either of the form of vertical accretions 
(Kleinecke, 1971; Sagar et al „ 1973, 1975) or change in the 
boundary conditions like river atage (Singh and Sagar, 1977). 

In this chapter, a numerical scheme for estimating 
aquifer parameters has been developed baaed on the analysis of 
historical record of aquifer response to the vertical 
accretions such as pumping, rainfall recharge and canal 
seepage. The principal permeability directions and coefficient 
of recharge due to rainfall can also, be estimated directly 
along with the aquifer parameters (T and S) using the inverse 
problem method. 

5.2 THE GOVERNING DIFFERENTIAL EQUATION 

The governing differential equation of two-dimensional 
transient ground water flow in a heterogeneous and isotropic 
confined aquifer can also be written as 

T + T *! h - + — ** -fA- ^ * Tyy ^ + Q = S 

XX £ X 2 yy <9x dx Sy Sy St 

(5.1) 

This equation is strictly valid for only confined 
aquifers as it assumes a time-invariant saturated flow and 
uniform horizontal velocity over the entire depth of the flow. 
However, this equation can be applied for unconfined aquifers 


I 
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If the following assumptions are satisfied. 

(1) Dupit-Forchhelmer ’ s assumptions namely 

(a) the velocity of the flow to be proportional to 
the tangent of the hydraulic gradient Instead of 
the sine as defined by Darcy’s equation, and 

(b) the flow to be horizontal and uniform everywhere 
in a vertical section (Todd, 1980). 

(2) Temporal fluctuations of water table, are small in- 
comparison to the mean saturated thickness. 


Sh Sh 
Sx 


) 


and 


Assuming the hydraulic gradients I > » 

ST ST ^ V 

transmissivity gradients J to be smal 1 , their 

products (Eqn.5.1) will be negligibly smaller in comparison to 

the other terms. Thus, neglecting the product terms, Eqn. 5.1 

reduces to the following form 


r *! h 

x * »x 2 


+ T 


* 2 h 




Q = S 


Sh 

St 


(5.2) 


vhere * - 2 % and are second spatial derivatives of the 
Sx Sy 2 

piezometric head in the directions of principal 
permeabilities. These can be derived In terms of their 
respective derivatives D p and D q in any two arbitrarily chosen 
orthogonal directions p and q as 


* h s D cos 2 0 + D sin £ 
* x 2 p q 


(5.3) 
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and 


© 2 h 


©y 


2 


D ain 2 © + D cos 2 © 
P q 


(5.4) 


In the above assumptions, the terms containing -s — -s — , 

Jp w 

Sh ae Sh se _ , ©h ©© 

©5“ ©x - ’ ©p~ 6y~ and ©q~©y have been neglected, because of 
their small order of magnitude. As shown in Fig. 5.1 © is the 
angle between two sets of orthogonal axes. Sagar (1975) 
presented this governing equation in terms of T , T and T 

pp qq pq 

2 

with the associated inequality constraint T + T -(T ) > 0. 

pp qq pq 

In the present approach, this constraint is eliminated and 
for the explicit estimation of the orientation of the 
principal permeability directions, equation 5.2 can be 
expressed in the form 

T f D cos 2 © + D sin 2 © 1 + f D sin 2 © + D cos 2 © 1 7+ Q = S 

xx t p q J L p q J ft \ at 

(5.5) 


For an anisotropic aquifer with spatial variation, the 
solution is obtained using the ‘direct approach’ (Neuman, 
1975) by solving this governing differential equation using 
finite element methods (Freeze and Cherry, 1979; Pinder and 
Gray, 1977). 

B. 3 DISCRETISATION OF THE GOVERNING DIFFERENTIAL EQUATION 

Equation 5.5 Is the continuous form of the governing 
differential equation. Its discretisation involves writing 
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down the following terms for discreet space and time points. 


i) Recharge (Q) 

ii) Spatial and temporal derivatives £- 
ili) Aquifer parameters (T and S) 


a 2 h a 2 h . ah 

~Z~2~ f - 2 and ~St 

Sq 


Employing the diacr et iaed form of Q a a obtained 

from Eqna .4.11 and 4.12, equation 5.5 can be expreaaed a a 

T xi( DP ik co * 2e + DQ ik oin2& ) + T yi ( DP ik * in2& + DQ ik coa2 ®] + 

f r [ P i,k’ P i,k-1’ P i,k-me’“] + 


C.. + X.. - U.. - S. DT.. = 0 

ik ik ik i ik 


(5-6) 


t b> 

where T x ^ , T yi and are aquifer parameters at i space 

point corresponding to T xx , T yy and S. And DP^* and 

tb 

are the discretised forms of DP, DQ and — ^ — for i space 
point at k time point. Eqns. 5.12 and 5.13 can be written as 


£ ( T xi’ T yi * & L’ S i ’ f r’ DP ik’ DQ ik' DT ik’ W ik’ X ik’ C ik’ 

P i,k ,P i ,k-l p i,k-me) = 0 


(5.7) 


Equations 5.6 and 5.7 express the governing differential 


equation in terms of 
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I) System properties 

T xi’ T yi ’ S i* ®i • 0,1 * f r 
il) Information derived from rainfall records 

P i,k’ P i,k-l’ ,P l,k-me 

lii) Information derived from piezometric head data 

( “i J , ( D ° ik ) . ( DT iJ 

iv) Information derived from hydrologic, agricultural 
and irrigation records 

( U ik] , ( C ik] , [ X ik) 

5. 4 . ESTI MATT ON OF DERIVATIVES 

The terms DP lk and DQ lk appearing In equations 5.6 and 
5.7 are estimates of the water table head In any two 
arbitrarily decided orthogonal directions p and q, at an angle 
& with principal permeability directions (p’and q’ ) as 
indicated in Fig. 5.1 and the term DT ik is the first temporal 
derivative. These derivatives can be directly obtained by 
differentiating the approxomating polynomial twice. Sagar et 
al . (1973, 1975) have proved that second spatial derivatives 

are necessary for estimating trensmlssivity by direct methods 
of Inverse problem. The least square polynomial apart from 
smoothening the data, also provides a closed form 

differentiable functional relationship between the water table 
head and the spatial coordinates. Thus, if H k (p,q) and 
H k are th * least a< l uare functional relations between 




Orthogonal axis (p) 

FIG. 5 1 DIRECTIONS OF PRINCIPAL PERMEABILITIES 



the water table head and the space coordinates p and q, at the 
beginning and end of the k^ time Interval, then 



C 5 . 8) 


DQ 


ik 




+ 




(5.9) 


DT 


H k+1 < p i ,q l* “ H k 


ik 


At. 


(5.10) 


where &re the coordlnetea (p , <1 ) o£ the i space 

point and At^ is the span of the k th interval. Equation 5.6 is 
the trlgonoalgebraic form of Eqn.5.2. The estimation of the 
derivatives (DP ik ), (DQ lk ), *nd (DT lfc ) can be made from 
historical water table data, and the estimates of ( u ik >> 
(C ), and (X,. ) from available historical data relating to 
th« groundwater withdrawals, canal seepage, geology, 
evapotran- splration etc.. Substitution of these estimates In 
Eqn . 5 . d for any given values of 1 and k. yields an equation In 
terns of syst.n paraneters. The equation thus evolved 1. 
indeterminate as It Involves (line) unknown terns. However, it 
can be converted to a deterministic one, when designed to 
estimate the properties of single space point which will 
consist of Eqn. 5. 6 written for 1 th space point and (dene) 
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different time periods i.e. f 


T xi[ DP ik co ‘ 2e 4 D « ik T yi [ DP lk eln 2 e 


+ DQ ik cosVJ + 


f r ( P i,k’ P i,k-1 P i,k-me’“ ] 


- S. DT.. = C.. + X .. - U.. 

i lk ik ik ik 


(5.11) 


where k = 1 (4+nc) and nc ia the number of 

rainfall-recharge parameters. 

The resulting system of non-linear simultaneous equations 
can be solved for (4+nc) unknown variables. This procedure can 
lead to many problems, since there are no fool-proof methods 
for the solution of non-linear simultaneous equations. The 
system of equations can be converted into a linear form by 
assuming linear rainfall relation and assigning some constant 
value to (from the known hydrogeological condition of the 
aquifer under study). 


S. 5 OPTIMIZATION 

Equation 5.6 is the discretised form of Eqn.5.2 and in 
addition, it incorporates the functional relationship for 
ralnfall-recharge. The right hand side of ' this equation may 
not be exactly equal to zero under the actual field 
conditions, due to any of the following reasons : 
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1) Ther* are many built-in assumptions in Eqn.5.2 which may 
not be always true. Neuman (1973) hae pointed out large 
number of physical situations which may violate many of 
these assumptions. These physical situations Include 
existence of fractures and dykes, unsaturated flow, 
compressibility, three- dimensional geometry, changes in 
the fluid density, non- Darcian flow in the areas of 
high Reynolds’ number. 

2) Numerical errors associated with the estimation of 


spatial and temporal 

derivatives 

of water 

table 

head 

CDP lk ),(D0 lk ) 

and (DT ik ) 

These 

errors 

can 

originate from the 

numerical 

algorithms 

adopted 

to 


arrive at the derivatives. 

3 ) The assumptions incorporated in the adopted rainfall- 
recharge relation (Eqns. 4.8, 4.9 and 4.10). 

4) The assumptions involved in calculating seepage from 
canals (the constants in the empirical formulae) may not 
exactly represent the study area. 

Because of these possible discrepancies, the solution of 


simultaneous equation may yield grossly misleading values of 
the parameters, and can not be used to estimate the parameters 
in most of the real world situations (Kashyap, 1981). 
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Keeping this In view, equation 5.6 can be expressed as 


T xi( DP ik col,2d + DQ ik sin2d ] + T yi [ DP ik sin 2 © + DQ ik cos 2 ©] + 

f r ( P i,k’ P i,k-1 P i,k-oe’°‘ 1 ] 


C ik + X ik ' U ik " S i DT ik 


E 


ik 


(5.12) 


where E^ k is the residue of the governing differential 
equation when written in the discrete form as given in 

Xi. x, C 

equation 5.6, for the 1 space point and k x time period. The 
exact magnitudes of (E^ k ) are not known. However, these 
residues will have certain frequency distribution. One of the 
most obvious frequency distribution for large sample is normal 
distribution with zero mean. The frequency distribution could 
be written as 


f (?) = - — 

/ 2 ncr 



(5.13) 


where o- is standard deviation and ? is frequency 
distribution. 


If this assumption relating to the normal distribution of 
the residues holds good, then the minimization of the sum of 
the squares of residues yield the most likelihood estimation 
i.e. minimize Y given by 
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y 


1 


n 

= E 

1=1 


m 

E 

k=l 



(5.14) 


The above distribution la valid only if sample size is 
large and there are no systematic errors. The other commonly 
used residue functional is the sum of moduli of the residues. 
The minimization of functional yields the moat likelihood 
solution if the residues are distributed as for the following 
distribution: 


f (?) = — 

Zncr 


exp 



(5.15) 


where is confidence level at p. 

Now the most likelihood estimates correspond to 


n m 

minimize E E l E itl 
1=1 k=l 


(5.16) 


The objective function to be minimized can Incorporate 

residue function by minimizing either squares of the residue 

or modulus of the residue, which are summed up over 

time-domain or over both time and space domains. Thus, the 

objective function can assume one of the following forms, 

apart from those given in Eqns.5.14 and 5.16. 

m „ 


2 

Y, 


E 

k=l 


m 


= E |E ik | 

k=l 


(5.17) 

for a given i, i = 1,2, . .n 
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where n la total number of apace points for which the 
parameters are to be estimated and m la the total number of 
time periods for which the data are available. 

a function of aquifer parameters and 
hydrogeological properties (i.e. the system parameters) of the 
I space point under consideration. The total number of 
decision variables Involved in the minimization of Y 2 equals 

the number of time periods i.e. k = 1 m. Thus, in each 

separate call to the optimization routine, aquifer parameters 
for one space point are estimated for different time periods 
(k). This approach allows number of decision variables to a 
manageable limit, but does not permit the Inclusion of 
constraints relating to the spatial variation of the aquifer 
parameters. This type of constraints may be necessary from the 
view point of getting a smooth solution (Emsellem and 
deMarsily, 1971). These can be Incorporated in the estimation 
of parameters provided the parameters of all the space points 
are estimated by the minimization of single objective 
function. This is permitted by minimization of (Eqns. 5.14 
and 5.15) which is function of aquifer parameters at all the 
space points. The total number of decision variables involved 
In the minimization goes up to (n X m) times for Y 1 . This 
increase in turn increases uncertainties generally associated 


with non-linear programming. 
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5.5.1 Constraint* 

The actual frequency distribution of the residues may 
follow the general pattern of the normal or exponential 
distribution, but may not strictly follow the geometrical 
relation as stipulated by the distributions (Eqns. 5.13 and 
5.15). As a result, the residues corresponding to the 
minimized residue functionals will seldom be identical to the 
actual residues. Thus, there may be an infinite number of 
near-optimal solutions which may vary considerably from the 
optimal solution but may lead to a much better predictive 
model (Neuman, 1973). In these cases, the knowledge relating 
to the geohydrology of the aquifer helps to choose the best 
set of parameters for obtaining the near-optimal solutions. 
These values of the parameters with judicious quantification 
can be incorporated in the parameter estimation program in the 
form of appropriate constraints which apart from ensuring 
non-negativity of parameters where required, may also 
stipulate the permissible range of variation or each of few of 
the constraints. The stipulated ranges of constraints may be 
more restrictive in nature, as the range is fixed on the basis 
of theoretical conditions. Similarly constraints can be 
imposed to restrict the spatial variation of. aquifer 
parameters and to obtain a ‘smooth solution' (Emsellem and 


deMarsily, 1971). 
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5.5.2 Sequential Unconstrained Minimization Technique 

The algorithm that can be employed to arrive at the 
optimal solution depends on the nature of the residue 
functional and constraints. The unconstrained minimization of 
the sum of squares of the residues can be accomplished by 
classical least square approach. This approach, widely used. 
In general involves a set of linear simultaneous equations. 
The constrained minimization of the sum of the squares of the 
residues or the sum of moduli of the residue can be carried 
out employing, with necessary modifications, one of the 
standard non-linear optimization algorithms. In the present 
work Sequential Unconstrained Minimization Technique, SUMT 
(McCormic and Flacco, 1968) has been used for the same. This 
program finds the minimum of a multivariable non-linear 
function subject to inequality and equality constraints. 

The SUMT uses the problem constraints and the original 
objective function which is minimized by any appropriate 
unconstrained multivariable technique. The steps involved in 
the algorithm are as follows : 

1. A modified objective function is formulated consisting 
of the original function and penalty function with the 
form where r is a positive constant. As the algorithm 
progresses, r is revaluated to form a monotonically 
decreasing sequence > r-^ > As r becomes small 

under suitable conditions, P approaches F and the 
problem is solved. 



112 

2. Select a starting point (feasible or nonfeasible) and 
an initial value of r. 

3. Determine the minimum of the modified objective 
function for the current value of r using SUMT or any 
appropriate technique. 

4. Estimate the optimal solution using extrapolation 
f romulae . 

5. Select a new value of r and repeat the procedure until 
convergence criteria are satisfied. 

This program consists of a main routine (to control 
subroutines that do special purpose of SUMT) and four user 
supplied subroutines. 

5.6 ESTIMATION OF AQUIFER PARAMETERS 

Aquifer parameters for the study area are estimated by 
considering the data of 20 discrete sequential periods of six 
months duration each. This was selected based on easy water 
budgeting nature in it i.e., when a water year is divided into 
two parts, each part has its own water-balancing nature and 
hence is a relatively closed system to operate. In the present 
study, the periods from December to May and from June to 
November are taken as pre-monsoon and post-monsoon periods 
respectively. Accordingly, the 66 months between June ’78 and 
November’ 83 have been divided into 11 slabs of six months 
each in the first five year set. The last six months period is 
taken as the overlap period (Table 5.1). 



Table 5.1 : Slabs for Paiameter Estimation (Set 1) 


I I J 


sl.no. period 


1 . 

Jun. ’78 

to 

Nov . 

’ 78 

2. 

Dec. ’78 

to 

May. 

. 79 

3. 

Jun. ’79 

to 

Nov . 

. 79 

4. 

Dec. ’79 

to 

May. 

’80 

5. 

Jun. '80 

to 

Nov . 

•80 

6 . 

Dec. *80 

to 

May. 

*81 

7. 

Jun. '81 

to 

Nov. 

•81 

8. 

Dec. ’81 

to 

May. 

•82 

9. 

Jun. ’82 

to 

Nov. 

’82 

10. 

Dec. ’82 

to 

May. 

*83 

11. 

Jun . ’ 83 

to 

Nov . 

, ’83 


slab no. 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


Overlap period 



The following data are considered for the modelling : 

a) The coefficients of the least square polynomial 
approximating for each of the periods which are calculated for 
all the sixty six months as explained in Chapter 3 and 
derivatives are calculated (as detailed in Section 5.4). These 
derivatives are stored as column matrices 

(**®ik} and (^ikj for ith «pace point. 

b) The rainfall figures for all the rain gauge stations, 
during each k*'* 1 period and stored as a column matrix £ P 
for all space points. The rainfall distribution graph is 

presented In Fig. 1.2. 

c) The effective ground water recharge and withdrawal 
values at each observation point (1), and for each of the 
periods, k, as [ cj ^ [ X ik J and [ « ik ] • 

These data were employed to generate the following column 
matrices for a given l** 1 space point. 

(“id, Kd, Kd. ( p ik). ( c ik), ( x ik) “ d ("id 

Each matrix has 11 elements corresponding to the 11 
periods. Marginal evapotranspiration values near the river and 
the canal banks have been incorporated into 
following initial values are adopted for 


( “id • The 

the decision 



variables . 
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S i = 0.2000, T xi = 150.00 m 2 /day, T yi = 150.00m 2 /day 
_ 0*00 , d ± =0.15 mm/month, a* = 0 . 8 5mro /month 

and a 3 = 2 Omni /month ... (5.18) 

The initial values of , T y ^ and are taken from the 
pump test data. The initial value of oi ^ is assigned assuming a 

a 

15% carry-over effect from the previous period. The decision 
regarding initial value of 6, a* and oi* were made depending on 

the basin characteristics. 

By the above setup, the residues oif the dlscretlsed 
Bousslnesq’s equation (5.11) for the i*** space point 
corresponding to all the eleven periods are completely 
described by the input data matrices and the initial values of 
the seven decision variables. For the estimation of optimal 
parameters, the following residue functional is adopted. 


Y 


i 


11 

E 


k=l 


E 


2 

ik 


Where Y^ 



(5.19) 


(5.20) 


The minimization of Y^ with respect 
decision variables is carried out subject to 


to these seven 
the following 


constraints . 
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00. 3C 

l > 

S i 

> 

00.03 

5000 

> 

T xi 

> 

40.00 

5000 

> 

T yi 

> 

40.00 

nfZ 

> 

a i 

> 

00.00 

00.20 

> 

i 

*i 

> 

00.00 

01.00 

> 

i 

“2 

> 

00.80 

80.00 

> 

i 

a 3 

> 

00.00 


Total number of constraints are fourteen. 

In the trial runs, aquifer parameters are estimated at 30 

space points by minimizing subject to the above 

constraints, using a non-linear optimization technique (SUMT). 

The results are far from the actual field values for all the 

space points, since maximum number of S values (20 out 30) are 

touching the lower limit (0.03) and values of T are closer to 

2 

the upper limit (5000 m /day). Even when the initial value of 
S is reduced to 0.00, all the estimates of S were still closer 
to the lower limit. At this stage where all initial values and 
constraints were rechecked and the entire data were thoroughly 
examined, it was observed that net inflows [c ik +x ik ] are far 
greater than net outflows £ U ik J , the excess of [ c ik +x ik ] ov « r 
[ U ik] ia cauaed due to ad ditional contribution to storage 

DT ik] and th ® net in£lows ( X ik +C ik) are directl y related to 
the estimation of higher values of T xi and T This problem 

can be solved by writing the equation in the form 



c ik * x lk 


’l* - ( - s i M ik ) * 

[ T xi ( DP ik coe2 ®j + D O ik sln^j J 4 
T yi( Dp ik » ln2 « l + »Q ik co. 2 e 1 J J 


(5.21) 


The excess values of * n t * ie optimal 
solution are assigned to the storage term at each space point. 
It may be indicated here that the estimation of the parameters 
is influenced by the Initial values assigned to the each of 
the space points. Hence, the aquifer parameters had to be 
estimated through trial and error by assigning initial values 
to all the 69 space points and checking the estimates. 

To solve this problem, a sequence of space points are 
estimated one after another by considering the nearest space 
point as a next point to be estimated from the first point. 
This scheme also did not give any improvement in the 
estimates. So another scheme was adopted, where, the estimated 
values at the first spatial point are given as initial values 


to the nearest space point. As a cross check, the estimated 
parameters at selected space points are compared with the 
actual values obtained through pump tests at the respective 


points. It has been observed that 
agreement within 10 percent 


both the values are in 
range. The estimated 


transmissivity values at the end of May 1983 are presented in 


Table 5.2. The storatlvity values for all the space points are 
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Table 5.2 Tranamlaalvity and Storage Eatimatea (Set 1) 


well 

No. 

transmleslvlty 

eq.m/day 

storage 

(MCM) 

1 

56.84 

30.39 

2 

171.65 

35.02 

3 

59.10 

66.26 

4 

128.19 

35.62 

5 

149.46 

143.42 

6 

169.63 

41.54 

7 

50.77 

33.27 

8 

161.74 

45.99 

9 

50.17 

34.17 

10 

50.03 

33.54 

11 

191.44 

24.72 

12 

161.46 

43.42 

13 

57.35 

37.21 

14 

153.04 

32.18 

15 

56.15 

58.47 

16 

52.88 

27.00 

17 

54.93 

35.89 

18 

51.94 

69.04 

19 

161.46 

21.26 

20 

159.22 

30.71 

21 

165.67 

55.55 

22 

163.70 

62.96 

23 

55.30 

63.87 

24 

50.42 

94.02 

25 

53.25 

47.97 

26 

50.17 

81.59 

27 

153.38 

10.40 

28 

156.82 

31.16 

29 

164.02 

53.44 

30 

53.23 

53.10 

31 

54.60 

6.81 

32 

50.59 

13.28 

33 

56.67 

22.14 

34 

52.31 

21.42 

35 

52.63 

57.39 

36 

195.77 

52.84 

37 

60.86 

86.86 

38 

58.94 

57.54 

39 

52.15 

33.98 

40 

50.46 

26.26 


Contd . . . . 
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Tabl € 5.2 Contd . . . 


well 

no . 

transmissivity 

sq.m/day 

storage 

(MCM) 

41 

50.29 

52.54 

42 

423.94 

13.94 

43 

428.69 

31.52 

44 

403.36 

44.20 

45 

71.46 

46.09 

46 

50.41 

31.51 

47 

51.41 

43.15 

48 

53.04 

65.65 

49 

53.87 

71.72 

50 

54.25 

91.78 

51 

210.35 

71.17 

52 

50.84 

36.24 

53 

54.25 

25.91 

54 

71.46 

36.17 

55 

220.85 

29.81 

56 

232.81 

55.04 

57 

51.95 

59.40 

58 

169.74 

39.47 

59 

58.26 

50.34 

60 

57.60 

69.09 

61 

51.03 

74.90 

62 

71.77 

63.10 

63 

51.01 

93.17 

64 

50.33 

67.27 

65 

169.74 

177.42 

66 

50.22 

72.60 

67 

54.57 

23.53 

68 

53.24 

40.07 

69 

54.72 

54.45 



estimated to be around 0.160. The angle between the principal 
permeability directions is estimated to be around 0.78 

radians . 

Storage values are estimated at all the 69 observation 
well points. For this purpose, one constant value for the 
thickness of the weathered and fractured zone is taken for 
each of the geological formations based on limited drill hole 
information. For granite and quartz arenites, this value was 
estimated as 35m from the surface while for limestones, this 
was fixed as 30m. The saturated thickness values (equal to the 
difference between the aquifer zone thickness and water table 
head) obtained at each of the observation points used for the 
estimation of storage. The estimates obtained for the two 
five-year block periods are indicated in Tables 5.2 and 5.5. 
The values vary from 10 to 140 MCM for granites which the same 
is in the range of 15 to 160 MCM for quartz arenites. The 
estimates for the limestone formations have worked out to be 
in the range 63 to 103 MCM. If information on the variation of 
thickness of weathered and fractured zones in the formations 
within the study area were available, the estimates of storage 
could have been more accurate. Further, it may be indicated 
here that although this storage is available around each well 
point, the extent to which the gorundwater exploitation can be 
carried out depends upon other management constraints. 

The coefficients of recharge due to rainfall were also 
estimated at each observation point along with the 



transmissivity and atorativity values. For an uniform value of 
0.15 given aa input, the model output has indicated a spatial 
variation between 0.13 and 0.18 (Table 5.3). 

The necessary software for the estimation of aquifer 
parameters in the present study was developed by the author. 
Algorithm for the same is prersented in Appendix B. 

S. 7 VALIDITY OF THE MODEL 

The procedure of parameter estimation by optimization 
essentially involves arriving at such estimates of the aquifer 
parameters which consistently minimize the residues of the 
Bousainesq’a equation, given by equation 5.12. These residues 
can be expressed in their simplest form as 


E * 




+ 0 



(5.22) 


The three terms on the right hand side of the equation 
5.23 represent the net excess inflow rate per unit area. 
Positive values of the second derivatives of the water table 

head indicate a concave surface and hence excess of inflows. 

S 2 h 

This is represented by — *■ > 0 and >©£, w ^ ere d «*nd ar ® 

angles between concave surface and its normal. Similarly 

negative second derivatives indicate a convex surface and the 

d^h 

resultant excess outflow represented by — ^ < 0 and 


A 



IZ. Z 


Table 5 


3 Estimated Rainfall Recharge Coefficient 


well coefficient 
no. 


well coefficient 
no . 


1 

0.1545 

2 

0.1566 

3 

0.1642 

4 

0.1478 

5 

0.1623 

6 

0.1385 

7 

0.1856 

8 

0.1783 

9 

0.1611 

10 

0.1677 

11 

0.1684 

12 

0.1874 

13 

0.1588 

14 

0.1500 

15 

0.1456 

16 

0.1626 

17 

0.1602 

18 

0.1683 

19 

0.1624 

20 

0.1545 

21 

0.1510 

22 

0.1764 

23 

0.1675 

24 

0.1620 

25 

0.1445 

26 

0.1415 

27 

0.1594 

28 

0.1528 

29 

0.1582 

30 

0.1485 

31 

0.1569 

32 

0.1475 

33 

0.1732 

34 

0.1580 

35 

0.1781 


36 

0.1324 

37 

0.1589 

38 

0.1652 

39 

0.1520 

40 

0.1753 

41 

0.1528 

42 

0.1577 

43 

0.1735 

44 

0.1712 

45 

0.1521 

46 

0.1687 

47 

0.1745 

48 

0.1824 

49 

0.1750 

50 

0.1648 

51 

0.1686 

52 

0.1618 

53 

0.1539 

54 

0.1605 

55 

0.1470 

56 

0.1547 

57 

0.1726 

58 

0.1663 

59 

0.1539 

60 

0.1824 

61 

0.1510 

62 

0.1574 

63 

0.1538 

64 

0.1825 

65 

0.1543 

66 

0.1652 

67 

0.1615 

68 

0.1707 

69 

0.1586 


z«ro second derivative indicates a plane surface with a 


couplets balance between the horizontal inflows and outflows 

and the expression for that becomes *■ = 0 and — S. . The 

A 2 1 ^ 

Sx 

terra Q represents net vertical accretion rate per unit area. 
The terra S ■ ^ represents the storage increase per unit area. 
Thus, Q denotes the net imbalance between the total inflows 
(horizontal and vertical) and the change in groundwater 
storage. As per the continuity equation, this imbalance should 
be zero for all the time periods. However, the residues may 
never vanish completely because of inadequate quantification. 

These residues are a function of T xx , T yy , S and other 
parameters representing the orientation of principal 
permeability directions and rainfall-recharge relations. Since 
these parameters are known or assumed to be time invariant, 
the objective is to arrive at such estimates of parameters 
for a specific space point, which consistently minimize as 
much as possible this imbalance at all the corresponding space 
points and for all the time periods. The minimization of the 
residue functional can be reviewed as multiple period water 
balance with time invariant parameters representing aquifer 


and hydrologic properties. 

To aetrtaln thie. a second s.t of data (Dec.iab.r 1982 to 
May 1988) pertaining to th. watar table heada, canal a.epage, 


rain fall “recharge and well draft etc. are con.ld.r.d a. Input 


values to th. aodel.The « aonth. ar. divided In to 11 slabs 


as tabulated in Table 5.4. 



Table 5.4 : Slabs for Parameter Estimation (Set 2) 

sl.no. period slab no. 


1 . 

Dec. 

’82 

to 

May. 

•83 

1 

2. 

Jun. 

*83 

to 

Nov. 

’83 

2 

3. 

Dec . 

’83 

to 

May. 

’84 

3 

4. 

Jun. 

*84 

to 

Nov. 

’84 

4 

5. 

Dec . 

’84 

to 

May. 

•85 

5 

6 . 

Jun. 

*85 

to 

Nov. 

’85 

6 

7. 

Dec . 

’85 

to 

May. 

'86 

7 

8. 

Jun. 

’86 

to 

Nov. 

•86 

8 

9. 

Dec . 

*86 

to 

May. 

•87 

9 

10. 

Jun. 

’87 

to 

Nov. 

•87 

10 

11. 

Dec . 

’87 

to 

May. 

•88 

Overlap period 
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As already Indicated the water table data for thia 66 
month period are also processed by least square polynomial 
approximation (Chapter 3) and the estimated derivatives are 
stored. The recharge components are also computed by the 
procedure explained in Chapter 4 and the generated values of 
recharge are also stored. Using all these Inputs, the aquifer 
parameters are estimated for the second set by taking the same 
initial values and constraints as for the first set to test 
the sensitivity of the model. The estimated transmissivity 
values for this period are presented in Table 5.5. As can be 
seen from the results, the difference between the 
transmissivity values for first five year period (1978-83) and 
second five year period (1982-88) is minimal. As expected 
there is no variation in storativity estimates. Unfortunately 
the T and S values as obtained form pump tests are few and 
hence cross checking was possible only to a very limited 


extent . 

The accuracy of the estimates was checked by simulating 


water table heads for the two five-year periods (1978 to 1983 
and 1983 to 1988) usin» direct problem (Peaceman and Rachford, 
1955). The procedure Involved in this direct problem is 
presented in Appendix C. These simulated value, of water table 
elevations are compared with collected observation well data. 

in each of the two five year blocks, results estimated 


for on. year duration each (June 1981 to May 1982 and December 
1983 to November 1984) are presented as typical examples. 
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T&bX« 5.5 TranamJ.ssIvi.ty and Storage Estimates (Set 2) 


well transmissivity storage 

no. sq.m/day (MCM) 


1 

56.56 

2 

173.12 

3 

59.96 

4 

127.21 

5 

150.88 

6 

170.54 

7 

51.45 

8 

163.19 

9 

51.08 

10 

50.86 

11 

192.21 

12 

162.43 

13 

58.19 

14 

154.28 

15 

55.94 

16 

52.92 

17 

54.72 

18 

53.22 

19 

162.43 

20 

160.52 

21 

162.03 

22 

161.76 

23 

56.11 

24 

50.97 

25 

52.70 

26 

51.47 

27 

158.41 

28 

156.27 

29 

166.89 

30 

53.72 

31 

54.65 

32 

51.08 

33 

57.24 

34 

52.58 

35 

51.13 

36 

194.16 

37 

60.29 

38 

60.55 

39 

52.59 

40 

51.22 


31.26 

36.30 

64.13 

35.23 

140.62 

43.70 

34.12 

42.47 
37.72 
34.33 
40.64 
43.42 

39.92 
32.60 
61.37 

28.44 

36.47 
69.91 
22.39 
30.05 
58.17 

65.35 
65.67 
96.34 

52.93 
83.58 

10.70 

31.36 

53.44 

54.71 
6.59 

14.13 

22.54 
21.85 
57.83 

58.55 

84.56 
60.99 
34.23 

26.57 


Contd 



Table 5.5 Contd 


well transmissivity storage 

no. sq.m/day (MCM) 


41 

50.63 

42 

424.28 

43 

430.85 

44 

408.68 

45 

72.60 

46 

51.86 

47 

50.99 

48 

55.00 

49 

54.16 

50 

54.79 

51 

215.09 

52 

50.88 

53 

53.65 

54 

71.15 

55 

221.35 

56 

238.10 

57 

51.86 

58 

162.72 

59 

58.95 

60 

58.34 

61 

52.09 

62 

72.11 

63 

52.17 

64 

50.23 

65 

160.23 

66 

50.69 

67 

53.45 

68 

54.29 

69 

54.68 


54.53 

15.17 

35.34 

42.14 

54.29 

31.96 

44.64 

70.21 

74.79 
92.31 
68.39 
38.84 
26.20 

37.95 
30.58 
57.83 

59.95 
38.23 

50.25 

63.72 

83.72 

63.79 
107.80 

' 73.62 
162.93 
73.47 

30.26 
42.75 
58.81 




Estimates for the period July 1981 (Table 5.6) indicated that 
the error varies between 0.001m and 9.724m. Out of 69 apace 
points in each month about 18 points have indicated an error 
percentage of 5% and above. For each of the space points the 
error is represented as positive or negative sign depending on 
whether the observed value is higher or lower than the 
estimates. Results for one month (September, 1984) presented 
in Table 5.7 also revealed a similar trend of variation 
between the observed and simulated head data. In this period 
also the number of nodal points where the error is above 5% is 
similar. The pattern of errors is identical in both the cases 
for any space point (Tables 5.6 and 5.7). 

The deviation of estimated values from the observed data 
for the groundwater heads can be due to several factors. The 
accuracy of these estimates to a certain extent can be 
improved by reducing the grid size as well as by reducing the 
time interval. This is true for the entire region for all the 
nodal points. However, the high values of errors observed at 
certain nodal points are predominantly due to the fact the 
effect of topographic, geological and land use factors are not 
completely taken into consideration in the recharge 
estimations. For example, the coefficient of rainfall recharge 
is chosen as 0.15 in the present study for the entire region, 
reflecting an overall average value. Due to non-availability 
of appropriate data, the spatial variation of this coefficient 
for different situations within the basin could not be 



Table 5. 


6 : Results of Back Calculation (Set 1) 
Period : July 1981 


well observed predicted error 

no. head(m) head(m) (*) 


1 

133.75 

2 

129.05 

3 

135.20 

4 

126.20 

5 

132.60 

6 

134.75 

7 

131.20 

8 

128.10 

9 

118.90 

10 

105.10 

11 

109.15 

12 

133.35 

13 

118.20 

14 

118.84 

15 

119.95 

16 

119.10 

17 

110.10 

18 

111.35 

19 

117.60 

20 

110.35 

21 

138.05 

22 

120.40 

23 

117.00 

24 

103.70 

25 

101.85 

26 

110.90 

27 

93.20 

28 

115.50 

29 

107.00 

30 

109.88 

31 

97.17 

32 

106.60 

33 

114.80 

34 

107.50 

35 

114.05 

36 

102.95 

37 

108.40 

38 

96.94 

39 

107.70 

40 

101.10 


131.08 

1.99 

127.48 

1.21 

134.84 

0.25 

129.27 

- 2.43 

138.01 

- 4.08 

127.02 

5.73 

128.19 

2.28 

123.93 

3.25 

115.54 

2.82 

109.78 

- 4.45 

107.18 

1.79 

134.74 

- 1.04 

124.66 

- 5.46 

120.07 

- 1.03 

122.34 

- 1.99 

118.99 

0.08 

113.05 

- 2.68 

108.94 

2.16 

109.41 

6.96 

118.32 

- 7.22 

128.79 

6.70 

125.62 

- 4.34 

117.06 

- 0.05 

110.63 

- 6.69 

103.94 

- 2.06 

103.24 

6.89 

100.93 

- 8.29 

113.42 

1.79 

107.53 

- 0.49 

110.92 

- 0.95 

109.21 

- 12.39 

103.51 

2.89 

108.51 

5.47 

104.38 

2.89 

112.69 

1.18 

106.85 

- 3.79 

101.57 

6.29 

98.04 

- 1.13 

98.52 

8.52 

104.38 

- 3.25 


. Contd. . 
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Table 5. 


7 : Results of Back Calculation (Set 2) 
Period : September 1984 

well observed predicted error 

no. head(m) head(m) (%) 


1 

134.25 

2 

130.65 

3 

134.55 

4 

125.60 

5 

131.50 

6 

134.50 

7 

131.85 

8 

131.50 

9 

120.90 

10 

106.10 

11 

110.85 

12 

132.70 

13 

119.62 

14 

121.94 

15 

120.90 

16 

121.10 

17 

110.15 

18 

110.35 

19 

120.70 

20 

111.15 

21 

137.85 

22 

120.15 

23 

117.17 

24 

104.95 

25 

103.85 

26 

111.80 

27 

95.70 

28 

114.90 

29 

105.80 

30 

110.65 

31 

101.90 

32 

109.05 

33 

115.15 

34 

109.10 

35 

114.40 

36 

104.10 

37 

107.15 

38 

98.50 

39 

107.30 

40 

102.20 


131.82 

1.80 

128.99 

1.26 

134.13 

0.30 

129.46 

-3.08 

136.75 

-3.99 

126.66 

5.82 

129.52 

1.76 

125.74 

4.37 

116.57 

3.57 

110.07 

-3.74 

107.98 

2.58 

133.97 

-0.96 

126.34 

-5.61 

121.83 

0.08 

124.23 

-2.75 

120.28 

0.67 

114.22 

-3.69 

110.10 

0.22 

113.29 

6.13 

119.22 

-7.26 

127.52 

7.48 

124.91 

-3.96 

118.58 

-1.20 

111.66 

-6.39 

105.06 

-1.17 

104.54 

6.48 

102.18 

-6.77 

113.54 

1.17 

108.08 

-2.16 

112.38 

-1.57 

110.70 

-8.63 

104.60 

4.07 

109.87 

4.58 

105.31 

3.46 

114.10 

0.26 

106.96 

-2.74 

102.13 

4.67 

99.03 

-0.54 

99.37 

7.42 

105.31 

-3.05 


Contd 



Table 5 


.7 Contd ... 

well observed 
no - head (m) 


41 

103.00 

42 

89.20 

43 

95.50 

44 

106.90 

45 

91.20 

46 

100.00 

47 

86.90 

48 

93.90 

49 

100.90 

50 

96.10 

51 

84.20 

52 

90.20 

53 

95.24 

54 

92.90 

55 

82.50 

56 

75.45 

57 

86.60 

58 

77.90 

59 

94.20 

60 

88.70 

61 

87.20 

62 

75.80 

63 

82.32 

64 

71.90 

65 

72.55 

66 

58.25 

67 

62.20 

68 

62.10 

69 

63.55 


predicted error 

head(m) (%) 


109.67 

- 6.47 

91.00 

- 2.02 

93.78 

1.79 

101.40 

5.14 

90.36 

0.91 

92.85 

7.14 

93.64 

- 7.76 

96.13 

- 2.37 

102.18 

- 1.26 

94.71 

1.44 

87.45 

- 3.86 

90.25 

- 0.06 

92.84 

2.51 

92.85 

0.04 

86.32 

- 4.63 

83.07 

- 10.11 

84.39 

2.54 

81.73 

- 4.92 

86.32 

8.36 

80.75 

8.96 

83.11 

4.68 

81.08 

- 6.97 

78.78 

4.29 

78.21 

- 8.77 

72.73 

- 0.25 

58.25 

- 0.00 

63.90 

- 2.74 

61.25 

1.36 

62.97 

0.90 



incorporated into the model. An enalysi. of date indicate, 
that the nodal pointe with maximum errora are altuated, in 
general , either cloeer to the canal areaa or in the vicinity 
of dykes. For example the nodal points 31, 37 and 5d are 

cloeer to the dyke, within the yranlte. Similarly the nodal 
points 27,60 and 67 are located close to lineaments observed 
on the satellite imageries. Nodal points 19, 20, 21, 27, 44, 

47 and 56 are close to the canal system in the study region. 
Although, for nodal points such as the ones described so far 
the effect of a dyke or a canal or a weak zone appears to 
influence the deviation of the simulated head values, the 
cumulative effect of various parameters including land slope 
and land use is responsible for the observed errors. 

To examine the influence of spatial variation of input 
data over the errors, split sampling technique has also been 
adopted. The coefficients of recharge and transmissivity 
estimates for various nodal points obtained as model ouput for 
the first five year block have been used as inputs for the 
respective nodal points for the simulation of water table head 
values for the second five year block. The deviation of the 
estimated head value from the observed head value at each of 
the observation wells has been calculated. The errors obtained 
through this approach for a period of one month (December 
1983), are compared with the errors worked out on the basis of 
back calculation for the same period (Table 5.8). As can be 
seen from the Table, there is a significant reduction in the 
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Table 5.8 : Water Table Head Values (m) Estimated by Split 

Sampling 


Period : December 1983 

well observed procedure A procedure B 

no. value estimated error(t) estimated error(%) 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 


133.56 
129.13 

133.41 

124.98 

130.15 
132.79 

130.53 

129.76 

119.36 

105.43 

109.27 
133.03 

118.21 

118.34 

120.37 

119.56 

109.54 
109.68 

118.25 
109.07 
137.09 
119.60 
116.63 

103.83 
100.79 

110.42 

94.25 

114.67 

105.81 

109.83 

98.20 
108.05 
113.65' 

109.20 

114.28 
103.00 

106.45 

97.20 
107.33 
100.49 


130.40 

127.46 

132.92 
128.33 
135.75 

125.92 

128.40 

124.16 

115.43 
108.95 
106.61 

134.40 
125.04 
120.48 

122.72 

118.97 

112.97 

108.78 

110.78 

117.38 

126.73 

124.25 

117.43 

110.35 

103.38 
102.75 
100.66 

112.2 

107.3 

111.25 

109.39 

103.20 
108.45 

103.85 
112.60 
106.06 

101.44 

97.77 

97.95 

103.85 


2.36 

1.28 

0.36 

- 2.68 

-4.31 

5.17 

1.63 

4.31 

3.28 

-3.34 

2.43 

-1.03 

-5.78 

-1.81 

-1.95 

0.49 

-3.14 

0.81 

6.30 

-7.62 

7.55 

-3.89 

-0.69 

-6.28 

-2.57 

6.94 

-6.81 

2.10 

-1.48 

-1.29 

-11.4 

4.48 

4.57 

4.89 

1.46 

-2.97 

4.70 

-0.58 

8.73 

-3.34 


131.21 

127.98 

134.21 
126.52 

129.16 

127.86 

127.36 

125.98 

120.00 

104.26 
107.23 

135.16 

123.77 

121.43 

121.98 

120.00 

111.54 

109.88 

114.21 
115.65 

129.89 

118.36 

116.17 

105.35 

99.98 

103.77 

96.73 

115.28 

106.43 

112.78 
103.77 

105.67 

114.35 

104.47 
113.74 
104.82 

103.55 

96.37 

104.63 

99.67 


1.75 

0.88 

-0.60 

-1.23 

0.75 

3.70 

2.42 

2.91 

-0.53 

1.10 

1.86 

-1.60 

-4.70 

-2.61 

-1.34 

-0.37 

-1.83 

-0.18 

3.41 

-6.04 

5.25 

1.02 

0.39 

-1.47 

0.79 

6.01 

-2.63 

-0.53 

-0.58 

- 2.68 

-5.67 

2.19 

-0.61 

4.32 

0.46 

-1.76 

2.71 

0.84 

2.51 

0.81 
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Table 5.8 Contd 


well observed procedure A (erocedure B 

no. Value estimated error(%) estimated error(t) 


41 

100.60 

109.20 

42 

88.40 

90.30 

43 

95.52 

93.08 

44 

104.05 

100.71 

45 

92.41 

90.05 

46 

99.87 

92.28 

47 

85.30 

92.80 

48 

94.55 

95.40 

49 

101.75 

101.83 

50 

95.90 

94.44 

51 

82.84 

85.77 

52 

90,01 

89.70 

53 

94.55 

92.26 

54 

91.66 

92.47 

55 

82.60 

85.66 

56 

75.66 

82.14 

57 

85.91 

83.46 

58 

76.51 

80.03 

59 

86.65 

84.8 

60 

86.65 

80.32 

61 

87.05 

83.61 

62 

75.10 

80.33 

63 

79.61 

77.23 

64 

74.50 

78.00 

65 

71.95 

70.75 

66 

57.76 

57.96 

67 

59.20 

62.77 

68 

60.33 

59.94 

69 

64.12 

61.64 


- 8.55 

103.22 

- 2.61 

- 2.16 

89.76 

- 1.54 

2.54 

94.51 

1.05 

3.20 

103.9 

0.13 

2.54 

93.65 

- 1.34 

7.59 

92.83 

7.04 

- 8.79 

88.27 

- 3.48 

- 0.90 

93.82 

0.76 

- 0.08 

102.63 

- 0.87 

1.51 

95.22 

0.70 

- 3.54 

83.70 

- 1.04 

0.33 

90.46 

- 0.50 

2.41 

93.60 

0.99 

- 0.88 

91.38 

0.29 

- 3.70 

83.74 

- 1.38 

- 8.57 

77.92 

- 2.99 

2.84 

86.92 

- 1.18 

- 4.60 

78.89 

- 3.11 

2.06 

85.43 

1.40 

7.30 

81.73 

5.66 

3.94 

88.17 

- 1.29 

- 6.97 

79.87 

- 6.35 

2.98 

78.74 

1.08 

- 4.70 

76.09 

- 2.14 

1.65 

72.45 

- 0.70 

- 0.34 

57.82 

- 0.10 

- 6.04 

58.63 

0.95 

0.63 

58.98 

2.22 

3.86 

62.51 

2.50 


( Procedure A : Back Calculation Method 
Procedure B: Split Sampling Method ) 


I 


error percentage when the split sampling method is adopted. 
Ulth this method, the maximum error is around 7% and seven 
values were having an error above 5%. In contrast to this, the 
maximum error was around 11% and 18 values were above 5% when 
back calculation method was adopted. 



CHAPTER 6 


INTEGRATED PICTURE 

The atudy area forms & part of the Nagarjuna Sagar left 
canal command area with three rivers (Musi river on west, 
Paleru river on east and Krishna river on south) and Nagarjuna 
Sagar left canal on the north as its boundaries. The basin 
covering about 830 sq.km is irrigated through the canals and 
groundwater for 9 months in a year. 

The geological sequence is of granites, quartz arenltes 
and limestones. The granites covering about 70% of the area 
are unclassified and belong to Archean age. Three small 
pockets of schists occur within the granitic terrain. In the 
western half of the granitic terrain, several dolerite dykes 
are present with a NE-SU trend. The granites are overlain by 
the quartz arenltes which in turn are succeeded by the younger 
limestones both belonging to the Kurnool Super Group. These 
two sedimentary formations occur in the southern part of the 
basin. The lineaments and the faults located in the area, in 
general, are in conformity with the dykes in their trend, 
possibly coinciding with the regional tectonic trend. 

For the groundwater evaluation, modelling techniques have 
been resorted to, based on the inverse problem as well as the 



direct problem. The former was adopted for the estimation of 
aquifer parameters while the latter was used for the water 
table head estimation to check the validity of the model and 
also for simulating the changes in water table elevation for 
particular well draft and canal seepages. 

The area has been divided into 20 x 19 grid (grid size 
2.5km x 2.5km). The water table heads pertaining to the study 
area were collected at 69 observation points and for 120 
months (June 1978 to May 1988). The data is divided into two 
blocks of 60 months each and each block is again divided into 
10 slabs of 6 months period (June to November, December to 
May). To remove the noise present in the water table data and 
to find out spatial and temporal variation of the ground water 
head, the data were processed by least square polynomial 
approximation. In this approximation, a functional approach 
was adopted by minimizing the sum of squares of residual 
errors. As a first attempt, the coefficients of third degree 
polynomial were calculated and the minimum standard error for 
a typical example was estimated as 0.9012m and variance value 
as 0.8910 after deletion of 15 space points and three 
polynomial terms. As a next step to avoid loss of data a 
fourth degree polynomial was considered and the coefficients 
of the same were calculated with 13 terms and 63 space points 
with a standard error of 0.7243 and a variance of 0.9714. 

It was observed that some space points in the region of 
dykes show anamalous values of groundwater head. An elevation 
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difference upto 15m has been evidenced for the water table 
heads at space points near dyke location. The polynomial 
coefficients generated by processing water table heads were 
used to calculate the second spatial derivatives and first 
temporal derivatives to help in aquifer parameter estimation. 
The software for processing water table data was developed by 
the author and the algorithm is presented in Appendix A. 

Water table contour maps prepared on the basis of 
observed head values reduced to mean sea level have clearly 
indicated the groundwater table configuration. The flow is 
from north with gradient in all directions. In general, an 
Increase in the spacing of the contours is evidenced in the 
vicinity of a dyke, lineament or a fault. The modification of 
the trend of contours within the study area is however due to 
a cumulative effect of all these structures as well as the 
canals which play an important role in their contribution to 
groundwater. The contour maps prepared indicating depth to 
water table from ground level have enabled a clear demarcation 
of water-logged localities in the study area. Trend surface 
maps also have been prepared to show the three-dimensional 

view of the groundwater table. 

Monthly rainfall data collected at two stations in the 
study area have been assigned to different observation points 
in the model for computation of rainfall recharge at each of 
the observation pooints. For this purpose, a value of 0.15 has 
been taken for the coefficient of recharge due to rainfall. 
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Canal seepage values are calculated as ratio of the 

difference between discharges at upstream and downstream sides 

of a canal to the area of its wetted surface. The recharge due 

to canal seepage was computed using the prescribed norms as 

outlined earlier (Chapter 4, Section 4.6) and the values are 

assigned to all the 69 observation points and for 120 months. 

The variation of recharge due to canal seepage at different 

3 3 

space points ranges for a typical month from 67x10 m to 
around 11738xl0 3 m 3 (Table 4.1). 

Effective withdrawals from all existing wells of the 
study area were estimated for all 120 months using the 
available norms and these values are also assigned to each of 
the observation points. 

The other components of groundwater balance equation such 
as evaporation and Irrigation effect were also considered. 
Monthly evaporation data were obtained and a single value is 
assigned to all the space points. 5% of the canal discharge 
values were added to the evaporation data towards irrigation 

effect. 

The inflows and outflows together with the net recharge 
values were calculated for all the months. Typical values for 
a month (March 1981) have been presented for all the 69 
observation points (Table 4.1). The values of recharge vary 

between 45 x 10 3 m 3 to around 11715 x 10 3 m 3 . 

Using all the values of inflows and outflows and 
derivatives of water table head, estimation of aquifer 
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parameters has been carried out by discretizing the 
Boussineseq’s equation and with the help of a non-linear 
optimization routine 'Sequential unconstrained Minimization 
Technique’, with chosen Initial values and constraints. The 
software for the same was developed by the author and the 
algorithm for the same is presented in Appendix B. Storativlty 
estimates were around 0.160, characteristic of the water table 
aquifers. Uhlle the transmissivity estimates in general were 
around 50 to 60 m 2 /day, locations with dykes, lineaments and 
faults have registered higher transmissivity values. In 
addition to the aquifer parameters, the angle between the two 
principal permeability directions and coefficient of recharge 
due to rainfall were also estimated through Inverse problem 
approach. The angle was estimated as 0.78 radians and the 
range of coefficient as 0.132 to 0.187. The storage values 
estimated at each of the observation wells ranged between 40 
MCM to 180 MCM. 

For checking the accuracy of the model, groundwater head 
values have been simulated using the estimated aquifer 


parameters 

and other 

related 

Inputs , 

by 

adopting 

the 

Alternating 

Direction 

Implicit 

method . 

The 

estimated 

head 


values are then compared with the actual field values for all 
the 120 months. It has been observed that the difference 
between these two sets of values varied from 0.5% to around 
11%. Of the 69 observation points, around 18 points were found 
to have an error above 5%. These anomalous values have been 
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attributed to the assumption in the calculation of recharge 
due to rainfall and to the effect of variation in geology, 
presence of dykes and land use. The errors in the estimation 
of water table heads could be further reduced significantly by 
assigning the transmissivity and the recharge coefficient 
estimates as obtained for all the 69 observation points at the 
end of the first five-year period, as inputs in the estimation 
of heads for the second block (1983-88). This split sampling 
procedure adopted had confirmed the effect of the spatial 
variation of the of the recharge due to geology and land use. 

APPLICATION OF THE MODEL FOR THE STUDY OF WATER LOGGING 

In the study area, since the inception of the canals the 
water table has registered an Increase. In three regions, the 
water level is up to 3m below the ground level (Fig. 6.1). In 
the region around Ponugodu in the NU part of the study area, 
the water table is about one metre from the ground resulting 
in water-logging conditions. To control this water-logging 
condition, a pilot study was conducted. This study involves 
introduction of a hypothetical drainage channel in the study 
area with assumed values of additional pumpage from the 
aquifer and simulating the aquifer behaviour by Alternating 
Direction Implict method (Peaceman and Rachford, 1955). This 
hypothetical channel has a length of about 16km and runs from 
observation point number 22 to observation point number 57 and 
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Joins the c&nal number 17 (Mutyala). The assumed discharge of 
this unlined canal is 2.832 cumecs. 

The data used for this study consists of the data at all 
space points comprising of the water table heads, storage 
factor, normal withdrawal rates and assumed additional draft 
values, transmissivity values and recharge factors etc. for 
the period December 1987 to May 88. The calculation of storage 
factor and value of error is similar as in the case of direct 
problem (Appendix C). The model is tested for five years and 
the head distribution at the end of May 1993 was estimated 
(Table 6.1). It is observed that with this assumed discharge 
values, the water table elevation in the study area would in 
general go down by about 1.0m at the end of five year period 
(May 1993). Particularly in the water-logged area (around 
observation point 21) in the vicinity of Ponugodu, decline in 
water table elevation would be as much as 3.00m. The contour 
plots for depth to water table for the month of December 1987 
are presented in Fig. 6.1 and the situation after introduction 
of the channel by May 1993 as simulated is shown in Fig. 6. 2. 

Thus, in the present work, groundwater modelling has 
enabled evaluation of groundwater stiuatlon in the canal 
command area in terms of transmissivity and storatlvity. 
Prediction of groundwater heads has also been made for 
particular well drafts and canal flows. This model can thus be 
used for monitoring water-logged areas and to arrive at 
solutions for the same. It is however, suggested that the 



144 


Table 6.1 Compariaion of Groundwater Levels (below ground 

surface) with Additional Pumpage 

well level in level in 

no. December’ 87 May’ 93 


1 

3.33 

3.56 

2 

4.65 

4.70 

3 

4.75 

5.70 

4 

3.80 

4.30 

5 

2.90 

3.11 

6 

4.75 

5.88 

7 

4.25 

5.65 

8 

5.30 

6.17 

9 

4.90 

5.07 

10 

3.50 

4.56 

11 

4.40 

4.85 

12 

3.70 

3.71 

13 

4.89 

6.02 

14 

11.93 

12.64 

15 

5.30 

6.88 

16 

4.50 

5.24 

17 

4.50 

5.85 

18 

4.50 

5.74 

19 

5.95 

6.61 

20 

5.02 

7.03 

21 

1.81 

4.63 

22 

5.28 

8.37 

23 

2.46 

5.55 

24 

4.22 

6.65 

25 

5.15 

6.72 

26 

6.42 

7.79 

27 

5.34 

6.92 

28 

3.55 

5.89 

29 

6.05 

8.82 

30 

4.20 

6.06 

31 

8.25 

9.11 

32 

6.65 

7.93 

33 

5.75 

6.01 

34 

4.55 

4.86 

35 

3.85 

4.31 

36 

6.03 

7.65 

37 

5.42 

7.86 

38 

7.00 

8.81 

39 

4.20 

5.64 

40 

3.70 

4.51 


Contd 



Table 6.1 : Contd. . . . 


Uell level in level in 

No. December’ 87 May’ 93 


41 

8.85 

9.65 

42 

5.35 

6.74 

43 

3.40 

4.10 

44 

4.94 

5.84 

45 

8.00 

10.28 

46 

4.20 

7.61 

47 

8.00 

9.25 

48 

4.75 

6.16 

49 

6.65 

7.52 

50 

8.10 

9.47 

51 

8.20 

9.29 

52 

4.40 

5.87 

53 

4.50 

5.61 

£4 

3.10 

3.49 

55 

4.57 

6.77 

56 

3.56 

4.24 

57 

3.38 

4.30 

58 

6.90 

4.12 

59 

12.70 

11.57 

60 

4.10 

5.05 

61 

10.50 

9.32 

62 

4.20 

4.67 

63 

8.60 

8.77 

64 

4.50 

5.06 

65 

8.40 

8.75 

66 

11.15 

11.46 

67 

18.25 

17.11 

68 

7.15 

7.71 

69 

4.25 

4.38 



1 *-4 w 


accuracy of the model can further be Improved by taking into 
account the geology, land use, topography and field features 
while assigning the rainfall recharge inputs to the model. In 
addition, it is also suggested that the measurement of water 
table heads is to be carried out, in a field problem like this 
type, with specific care. It has to be ensured that when the 
water level in a well is measured, there is no interference of 
the same due to pumping at the adjacent sites. 








REFERENCES 


Abouflrassl, M. , and Marino, M.A. , 1984. Cokriging of 
aquifer transmissivities from field measurements of 
transmissivity and specific capacity, Math.Geol., 16(1) 

A. P -S . G. U. D . REPORT, 1981. Quarterly report on 
groundwater situation in Nalgonda district, Andhra. Pradesh. 
State. Ground. Uater .Department., Hyderabad. 

Bard, Y. , 1974. Nonlinear Parameter Estimation, John 
Uiley , New York. 

Bear, Jacob, 1972. Dynamics of Fluids in Porous Media, 
American Elsevier Publ. Co. Inc., New York. 

Beck, V. J, and Kenneth, A. J. ,1977. Parameter Estimation 
in Engineering and Science, Uiley Interscience, New York. 

Bedinger, M.S. , Red, J.E. and Griffin, J.D., 1973. 
Digital Computer Programs for Analysis of Groundwater flow, 
U.S.Geol. Survey open file report. 

Bellman, R. and Kabala, R., 1965. Quasilinearization and 
Non-linear Boundary Value Problems, American Elsevier 
Publ. Co., Inc., New York. 

Becker, L., and Yeh, U.U.G., 1972. Identification of 
parameters in unsteady open-channel flows, Uater Resour. Res., 
8(4) , 956-965. 

Blrtles, A.B., and Morel, E.H.,1979. Calculation of 
aquifer parameters from sparse data, Uater Resour. 
Res. ,15(4), 832-844. 


Bittenger, M.U. , Duke, H.R. and Longenbaugh, R.A., 1967. 
Mathematical simulation for better aquifer management, 
Publication 72, IASH Symposium, Haifa, 1967, 509 519. 


Bredehoeft, J.D. and Pinder, G.F., 1970. Digital 
of areal flow in multi aquifer groundwater systems : 
three dimensional model, Uater Resour. Res . , 6 ( ), 


analysis 
a quasi 
883-888. 


Britles, A.B. 
aquifer parameters 
15(4), 832-844. 


and Morel, E.H., 1979. Calculation of 

from sparse data, Uater. Resour. Res., 


Chang, S., and Yeh, U.U.G., 1976. 
for the solution of the large-seal 
groundwater, Uater Resour. Res., 12(3), 


A proposed algorithm 
e Inverse problem in 
365-374. 



Chaturvedi, R.S. and Chandra, S. , 1961. Analytical 
approach towards solving water logging problems. Research 
Journal, Vol . IV, No. 1, University of Roorkee, Nov. 1961. 

Chatter jee, S. , 1977. Regression Analysis by Example John 
Uiley and Sons, New YorK. 

Chavent, G. , 1974. Identification of functional 
parameters in partial differential equations, in 
Identification of Parameters in Distributed Systems, edited by 
Goodson, R.E. and M. Polls, pp. 31-48, American Society of 
Mechanical Engineers, New York. 

Chavent, G., 1979. About the stability of the optimal 
control solution of Inverse problems, in Inverse and 
Improperly Posed Problems in Differential Equations, edited by 
G. Anger, pp. 45-48, Akademi e-Verlag , Berlin. 

Chavent, G. , 1979. Identification of distributed 
parameter system: About the output least square method, its 
implementation, and ident i f iabi 1 ity , in Identification and 
System Parameter Estimation, edited by R. Isermann, vol.l,pp. 
85-97, Pergamon, New York. 

Chavent, G. , 1983. Local stability of the output least 
square parameter estimation technique. Math. Appl .Comp. , 2(1), 
3-22. 

Chavent. G. , Dupuy, M. and Lemonnier, P. 1975. History 
matching by use of optimal control theory, Soc . Pet . Eng . J . , 
15(1), 74-86. 

Chen, U.H., Gavals Gr. S., and Uasserman, M.L., 1974, 
Anew algorithm for automatic history matching, Soc. Pet. Engg. 
J., 14(6), 593-608. 

Coates, K.H., Dampsey , J . R . , and Henderson, J.H., 
1970. A new technique for determining reservoir description 
from field performance data, Soc . Pet . Eng . J . , 10(1), 66-74. 

Cooley, R. L., 1977. A method of estimating parameters 
and assessing reliability for models of steady state 
groundwater flow. 1. Theory and numerical properties, Uater 
Resour. Res., 13(2), 318-324. 

Cooley, R. L., 1979. A method of estimating 
parameters and assessing reliability for models of steady 
state groundwater flow. 2. Application of statistical 
analysis, Uater Resour. Res., 15(3), 603-617. 

Cooley, R.L., 1982. Incorporation of prior information on 
parameters into nonlinear regression groundwater flow models, 

1. Theory, Uater Resour. Res., 18(4), 965-976. 

Cooley, R.L., 1983. Incorporation of prior information on 
parameter into nonlinear regression groundwater flow models, 

2. Applications, Uater Resour Res., 19(3), 662-676. 



I U I 


Cooley, R.L., and Sinclair, P.J., 1976. Uniqueness of a 
model of steady-state groundwater flow, J.Hydrol., Vol.31, 

245 - 269 . 

Cooper, H.H., 1966. The equation of groundwater flow in 
fixed and deforming coordinates, Journal Geophysical 
Research, 71(20), 4785-4790. 

Crank, J. and Nicolson, P., 1947. A practical method for 
numerical evaluation of solutions of partial differential 
equations of the heat conduction type. Proc. Camb. Phil. Soc. 
Vol . 43, 50-67. 

Dagan, G. , 1985. Stochastic modelling of groundwater flow 
by unconditional and conditional probabilities: The inverse 
problem, Uater Resour. Res., 21(1), 65-72. 

Daniel, C. and Uood, I.S., 1971. Fitting Equations 
to Data, Uiley Interscience, 

Davis, J.C., 1973. Statistics and Data Analysis in 
Geology. John Uiley & sons New York. 

Delhomme, J.P., 1978. Kriging in Hydrosciences, 
Advances in Uater Resources, Vol. 1. No. 5. 

Delhome, J.P., 1979. Spatial variability and uncertainlty 
in groundwater flow parameters: A geostat 1st leal approach, 
Uater Resour. Res., 15(2), 269-280. 

Dewiest, R.J.M., 1966. On the storage coefficient and the 
equations of groundwater flow. Journal Geophysical Research, 
71(4). 


DIStefano, N. , and Rath, A., 1975. An identification 
approach to subsurface hydrological systems, Uater 
Resour. Res., 11(6), 1005-1012. 


Dogru, A.H., and Seinfeld, J.H., 1981. Comparison of 
sensitivity coefficient calculation methods in automatic 
history matching, Soc . Pet .Eng. J - , 21(5), 551-557. 


Dogru, A.H., Dixon, T.N., and 
Confidence limits on the parameters 
slightly compressible, single 

Soc .Pet . Eng. J - , 17(1), 42-56. 


Edgar, T.F., 1977. 

and predictions of 
phas e res er vo i rs , 


Dooge, J.C.I., 1973. Linear theory 

hydrological systems. Agricultural Research 

United States Department of Agriculture, 

Bulletin No. 1468. 


of 

Service. 

Technical 


Douglas, J. Jr. , 1962. Alternating direction methods for 
free space variables, Num. Math., Vol. 4, 41-63. 


Dreizln, 
of response 


Y.C. and Haimes, Y.Y., 1977. A 

functions for groundwater management 


hierarchy 

Uater 



Resour. Res. 13,(1), 78-86. 


Drury, S.A. and Holt, R.U. , 1980. The tectonic frame 
work of the South Indian craton a reconnaissance involving 
LANDSAT Imagery, Tectnophysics , Vol.65., T1-T15. 

Eduardo, A., Remson, I., Pikul , M.F. and Thomas, 
U.A. , 1974. Optimal pumping for aquifer dewatering. 
Journal of Hydraulic Division. ASCE, 100(HY7), 869-877. 

Emsellem, Y. and deMarslly, G. , 1971. An automatic 
solution for inverse problem, Uater Resour. Res. 7(5), 
1264-1283. 


Ferris, J.G., 1952. Cyclic fluctuations of water levels 
as a basis for determining aquifer transmissibillty , 
U.S. Geological Survey, Groundwater, Note 1. 


France, P.U., 1974. Finite element analysis of 
three dimensional groundwater flow problems, J. Hydrol . , Vol . 
21, 381-398. 


Freeze, R. A., and Cherry, J.A., 1979, Groundwater, 
Prentice-Hall Inc. 

Frlnd , E.O. and Plnder, G.F., 1973. Galerkln solution 
for the inverse problem for aquifer transmissivity, Uater. 
Resour. Res., 9(5), 1397-1410. 


Gates, J.S. and Kislel, C.C., 1974. Uorth of additional 
data to a digital computer model of groundwater basin, Uater 
Resour. Res., 10(5), 1031-1038. 

Gavalas, G.R., Shah, P.C. and Seinfeld, J.H., 1 J 76 « 
Reservoir history matching by Bayesian estimation, Soc . Pet. 
Eng., J-, 16(6), 337-350. 


Guitin, M.E., 1964. Variational principles for lin *J 

initial value problems, Quarterly Applied Mechanics, Vol. a 


2522 pp. 

Haimes, Y.Y. and Dreizln, 
groundwater and surface water 
Resour. Res., 13(1), 69-77. 


Y.C., 1977. Management of 

through decomposition, Uater 


Haimes, Y.Y., Perrlne, R.L. and Uismer, D.A. 1968 

Identification of aquifer parameters Resources 

multilevel optimisation, Contribution No. 123, Uater Resources 

Center., Ulnv. Calif., Berkeley. 


Hall, U.A. 
Resources Systems 


and Dracup, J.A., 
Engineering, McGraw-Hill, 


1970. 
New York. 


Uater 


Hefez, E . , Shamir, V. and Bear, 
the parameters of an aquifer cell model 
Res., 11(6), 993-1004. 


1975. Identifying 
Uater Resour. 



IDO 


Hoeksema, R.J., and Kitanidia, P.K., 1984. An application 
or the geostatlstical approach to the inverae ‘problem in 
two-dimenaional groundwater modeling, Water Reaour. Rea. 
20 ( 7 ), 1003-1020. 

Hoekaema, R.J., and Kitanidia, P.K., 1985. Comparialon of 
Gauaaian conditional mean and kriglng eatimation in the 
geostatlstical aolution of the Inverae problem, Water Reaour 
Rea., 21(6), 825-836. 

Jacob, C.E., 1950. Flow of groundwater. Engineering 
Hydraullca, H.Rouae (Ed.), New York, John Wiley and Sona, 
New York. 

Jacquard, P., and Jain, C. , 1965. Permeability 

diatrlbution from field preaaure data, Soc. Pet. Eng. J., 
5 ( 4 ), 281-294. 

Jahns, H.O., 1966. A rapid method for obtaining a 
two-dimenaional reaervolr deacrlptlon from well preaaure 
response data, Soc. Pet. Eng. J. , 6(4), 315-327. 

Javendel , I. and Witherspoon, P.A., 1968. Analysis 

of transient fluid flows in multi-layered systems,. 
Contribution No. 124, Water Resources Centre, Univ. Calif., 
Berkel ey. 

Kafri , U. and Aaher , J.Ben, 1978. Computer estimates of 
natural recharge through soils in southern Arizona, U.S.A. 
J . Hydrol . , Vol . 38, 125-138. 

Karanth, K.R., 1987. Groundwater Assessment 

Development and Management, Tata McGraw Hill Publishing Co. 
Ltd., New Delhi. 576-657. 

Kaahyap, D. , 1981. Mathematical Modeling of groundwater 
systems, Unpublished Thesis, School of Hydrology, Univ. 
Roorkee, Roorkee. 

Kaahyap, D. , and Chandra, S. , 1982. A nonlinear 

optimization method for aquifer parameter estimation, J. 
Hydrol., 57, 163-173. 

Kitamura, S. , and Nakagiri , S., 1977. Identifying of 

spatially-varying and constant parameters in distributed 
systems of parabolic type, SIAM J. Contr. Optimiz., 15(5), 
785-802. 

Kitanidia, P.K., and Vomvoris, E.G.,1983. A 

geostatist leal approach to the inverse problem in groundwater 
modeling (steady state) and one-dimensional simulations, water 
Resour. Rea., 19(3), 677-690. 

Kleincke, D., 1971. Use of linear programming for 

estimating geohydrologlcal parameters of groundwater basins, 
Uater Resour. Res., 7(2), 367-374. 



Kleincke, D. , 1972. Comments on "An automatic solution 
the inverse problem" by Y. Emsellem and G. deMarslly, 
Uater Resour. Res., 8(4), 1128-1129. 

Knowles, T.R. , Calborn, B.J. and Uells, D.M. , 1972. 
A computerised procedure to determine aquifer 
characteristics. Department of Civil Engineering, Texas Tech. 
Univ., Lubbock, Texas, URC-72-5. 

Krishnamurti , N. , 1977. Mathematical modelling of 
natural groundwater recharge, Uater Resour. Res. 13(4), 

720 - 724 . 


Kruseman, G.P. and De Rldder, N.A. , 1970. Analysis and 
evaluation of test pumping data. International Institute for 
land reclamation and improvement. Uageninger, The Netherlands . 

Kubrusly , C. S., 1977. Distributed parameter system 
identification, a survey, Int. J. Contr. , 26(4), 509-535. 


Labadie, J.U., 1972. Decomposition of large scale 
nonconvex parameter identification problem in geohydrology. 
Rep. ORC 72-73, Oper. Res. Cent., Univ of Calif., Berkeley. 


Labadie, J.U., 1975. A surrogate parameter approach 
to modelling groundwater basins, Uater Resources 
Bulletin, 11(1), 97-113. 


Lakshminarayana, V. and Rajagopalan, S.P., 1978. Digital 
model studies of unsteady state radial flow to partially 
penetrating wells in unconfined anisotropic aquifers, J. of 
Hydrol . , Vo 1 - 38, 249-262 . 

Lin, A.C., and Yeh, U.U.G., 1974. Identification of 
parameters in an Inhomogeneous aquifer by use of the maximum 
principle of optimal control and quaslllnearlzatlon, Uater 
Resour. Res., 10(4), 829-838. 


Lovel , R.E., Duckstein, L. and Kisiel, C.C., 1972. Use 
of subjective Information in estimation of aquifer 
paramet ers , Uat er Resour. Res., 8(3), 680 690. 

Maddock, T., Ill, 1972. Algebraic technological function 
from a simulation model, Uater Resour., Res., 8(1) 129-134. 


Maddock T III 197 6. A drawdown prediction model based 
on regression L^lysli, Uater Resour. Res., 12(4) 818-822. 


Maddock, T. , III end Haimes, Y.Y., 1975. 
the planning and management of groundwater, 
Res. 11(1), 7-14. 


A tax system for 
Uater Reaourcea 


Marino, M. A. and Yeh, U.U.G., 
parameters in finite leaky aquifer 
Div. ASCE 99 (HY2 ) , 319-336. 


1973. Identification of 
system. Journal Hydr. 


for 



*«ri«?Ji qUar ? t ’ U ’’ 1963 * An algorithm 
estimation of nonlinear parameters, SIAM, 


for least squares 
J-, 11, 431-441. 


McCormick, G. and Fiacco, A., 1968. 

Programming Sequential Unconstrained Minimisation 

York. ’ 


Non-linear 
Ulley, New 


Monech A.F. and Kisiel, C.C., 1970. Application 

of convolution relation to estimating recharge from an 
ephemeral stream, Uater Resour. Res., 6(4), 1087-1094. 

Morel-Seytoux, H.J., 1975 . A combined model of water- 
table and river stage evolution, Uater Resour. Res., 11(6). 

Mor el — Seyt oux , H.J., 1976. Derivation of equat ion for 
rainfall infiltration, J.Hydrol., Vol.31, 203-219. 

Murray, U.A. and Johnson, R.L., 1977. Modelling of 

unconfined groundwater systems. Groundwater, 15(4), 306-312. 

Navarro, A., 1977. A modified optimization method of 
estimating aquifer parameters, Uater Resour. Res., 13(6), 
935-939. 

Nelson, R.U., 1962. In-place measurement of 

permeability in heterogeneous media, 1, Theory of proposed 
method, J. Geophys. Res., 65(6), 1753-1758. 

Nelson, R.U., 1968. In-place determination of 

permeability distribution for heterogeneous porous media 
through analysis of energy dissipation, Soc.,Pet., 33-42. 

Neuman, S.P. and Uitherspoon, P.A., 1970. Variational 

principles for confined and unconfined flow of groundwater, 
Uater Resour. Res., 6(5), 1376-1382. 

Neuman, S.P., 1973. Calibration of distributed 
parameters in goundwater flow models viewed as multiple 
objective decision process under uncertainty, Uater Resour. 
Res. 9(4)., 1006-1021. 

Neuman, S.P., 1975. Role of subjective value judgement in 
parameter identification. In G.C. Van Steenkiste, (ed.). 
Modelling and simulation of water resources systems, North 
Holland, Amsterdam. 

Neuman, S.P., and Yakowltz, 1979. A statistical 
approach to the inverse problem of aquifer hydrology, 1, 
Theory, Uater Resour. Res., 15(4), 845-860. 

Neuman, S.P., 1980. A statistical approach to the inverse 
problem of aquifer hydrology, 3, Improved solution method and 
added perspective, Uater Resour. Res., 16(2), 331-346. 

Nutbrown, D.A. , 1975. Identification of parameters in 
linear equation of groundwater flow, Uater Resour. Res., 11(4) 



Placeman, D.W. and Rachford. H.H., 1955. The 
numerical solution of parabolic and elliptic differential 
equations, SIAM.J., Vol.3, 28-41. 

Pinder, G.F. 1970. A Galerkin infinite element 
simulation of groundwater contamination on Long Island, New 
York, Water Resour. Res., 9(6), 1657-1669. 

Pinder, G.F. and Bredehoeft, J.D., 1968. Application of 
the digital computer for aquifer evaluation, Water Resour. 
Rea., 4(5), 1069-1109. 


Pinder. G.F. and Frlnd, E.O., 1972. Application of 
Galerkln's procedure to aquifer analysis. Water Resour. Res., 
8(1), : 108-120. 

Pinder, G.F. and Gray, W.G.,1977. Finite Element 

Simulation in Surface and Subsurface Hydrology, Academic 

Press, Inc., London. 


Pinder. G.G. , Bredehoeft, J.D. and Cooper, H.H. Jr., 
1969. Determination of aquifer diffusivity from aquifer 
response to fluctuation in river stage. Water Resour. 
Res., 5(4), 850-855. 


Prickett, T . A. and Lonnquist, C.G., 1971. Selected 
Digital Computer Techniques for Groundwater Resource 
Evaluation. Ill, State Water Survey Bulletin, Vol.55. 


Ralston, Anthony, 1965. A First Course in Numerical 
Analysis, McGraw-Hill Inc. 

Reeves, M. , 1975. Estimating infiltration for 

erratic rainfall. Water Resour. Res., 11(1), 


Remson, Irwin, Hornberger, G.M. and Molz, F.J., 1971 

Numerical Methods in Subsurface Hydrology, Wiley Interscience 


Richardson, L.F., 1910. The approximate arithmetical 
solution by finite differences of physical problems involving 
differential equations with an application to the stresses In 
a masonry dam, Phil. Trans. Royal Soc . , A 210, 307-357. (As 
referred In Kashyap , 1981) 


Rowe P.P., 1960. An equation for estimating 

transmlsslbi llty and coefficient of storage from J- iver 
fluctuations. Journal of Geophysical Research, 65(1U), 
3419-3424. 


Rushton, K.R., 1974. Critical analysis 

alternate direction implicit method of aquifer analysis, 
Hydrol . , Vol.21, 153-172. 


of 

J. 


Rushton, K.R. and Ward, C., 1979 

groundwater recharge, J. Hydrol., Vol-41, 

Rushton, K.R. and Rathod, K.S. , 


The estimation of 
345-361. 

1985. Horizontal 



and vertical components of flow deduced 
heads, J.hydrol., Vol.79, 261-278. 


from 


groundwat er 


Rushton , K.R., 1986. Groundwater Models in Groundwater: 
Occurrence, Development and Protection, Chapter 6. 
Ed. Brandon. T.U., Uater Practice Manuals, The Institution of 
Engineers and Scientists, London. 

Sadeghipour , J. and Yeh, U.U.G., 1984. Parameter 
identification of groundwater aquifer models: A generalized 
least squares approach, Uater Resour. Res., 20(7), 971-979. 

Sagar, B. and Singh, S.R., 1979. Aquifer diffusivlty from 
noisy boundary data. Journal of Hydr . Div . ASCE 
105(HY8), 943-954. 


Sagar, B., Kisiel, C.C. and Ducksteln, L. , 1973. 
Estimation of aquifer properties from field data. Proceedings 
International Symposium on Development of Groundwater 
Resources, Madras, India. No.l, Vol.I., 26-29 

Sagar, B. , Yacowitz, S. and Duckstein, L. , 1975. A direct 
method for Identification of parameters of dynamic 
non-homogeneous aquifers, Uater Resour. Res., 11(4), 563-570. 


Seethapathl, P.V., 1984. Tyson-Ueber Groundwater Flow 
Model, UM-I, National Institute of Hydrology, Roorkee. 

Shah, P.C. .Gavalas, G.R. and Seinfeld, J.H.,1978. Error 
analysis in history matching: The optimum level of 
parameterization, Soc . Pet. Eng. J., 18(3), 219-228. 

Singh, B.P. and Chandra, S. , 1977. Evaluation of the 
optimal thickness of soil between source and detector in the 
gamma ray transmission method. J. Hydrol . , Vol.32, 

189 - 191 . 

Singh, S.R. and Sagar, B., 1977. Estimation of aquifer 
diffusivlty in stream-aquifer systems, Journal of 
Hydr . Div. ASCE, 103(HY11), 1293-1302. 


Sokolov, A. A. and Chapman, T.G. (Eds.), 1974. 
Methods of Uater Balance Computations -An International Guide 
for Research and Practice, The UNESCO Press. 


Stallman, R.U. , 1956. Numerical analysis of regional 
water levels to define aquifer hydrology, Trans. AGU, 37(4), 
451-460. 


Stone, H.K., 1968. Iterative solution of Implicit 
approximation of multidimensional partial differential 
equations, Society of Industrial Applied Mathematics. Anal. 
Vol.5, No. 3, 530-558. 

Sukhija, B.S. and Shah, C.R., 1976. Conformity of 
ground water recharge by tritium method and mathematical 
modelling, J. Hydrol., Vol.30, 167-178. 



Sun, N.Z., and Yeh, W.W.G. 1985 
parameter structure In groundwater 
Reaour . Rea., 21(d), 869-883. 


Identification of 
Inverse problem, Water 


Theis, C.V., 1935. The relation between lowering 

of piezometric surface and the rate and duration of a well 
using ground water storage, Transactions Am. Geophy. Union 
Vol.16, 519-524. (As referred in Yeh, 1986). 

Thomas, L.K., Heliums, L.J., and Reheis, G.M. 1972. A 
nonlinear automatic history matching technique for reservoir 
simulation models, Soc.Pet.Eng. , J., 12(6), 508-514. 


TIson, G. Jr., 1965. Fluctuations of groundwater levels, 
Advances In Geophysics Vol.2, 303-325. 


Todd, D.K., 1980. Groundwater Hydrology, John Wiley and 

Sons, New York. 


Tomilson, L.M. and Rushton, K.R., 1975. The alternating 
direction explicit method for analysing groundwater flow, 
J. Hydrol . , Vol.27, 267-274. 

Trescott, P.C., 1975. Documentation of 
finite-difference three dimensional groundwater flow, U.S. 
Geol . Survey Open file Report 75. Trescott, P.C. and 
Larson, S.P., 1977. Comparison of iterative methods of 
solving two-dimensional groundwater flow. Water Resour. Res., 
13(1), 125-136. 

Trescott, P.C., Pinder, G.F. and Larson, S.P., 1976. 
Finite difference model for aquifer simulation In two 
dimensions, with results of numerical experiments , U.S. Geol. 
Survey, Techniques Water Resources Inv. , Book 8, Chap. Cl. 


Tyson, H.N. and Webor, E.M., 1964. Groundwater 
management for national future-Computer simulation of 
groundwater basins, Journal of Hydraulic Div. Proc. ASCE, 
Vol . 90 , No. 4, 59-77. 

Vemuri, V. and Karplus, W.J., 1969. Identification 
of non-linear parameters of groundwater basins by 
hybrid computations. Water Resour. Res., 5(1), 172-185. 


Venetls , C., 1971. Estimating infiltration and/or the 
parameters of unconflned aquifers from groundwater level 
observations, J. .Hydrol., Vol. 12., 161-169. 


Vogel, J.C., Thilo, L. 
Determination of groundwater 
J. Hydrol., Vol. 23, 131-140. 


and Van Dijken, M. , 1974. 

recharge with tirltlum, 


Walton, W.C., 1970. Ground Water Resources Evaluation, 
McGraw-Hill, Singapore . 



I u >J 


Uang, F.H. and Anderson M.P., 1982. Introduction to 
Groundwater Modelling, Finite Difference and Finite Element 
Methods , U.H. Freeman and Co. , Sanf ransisco . 

Uatson, P, Sinclair, P. and Uaggoner, R.K, 1976. 
Quantitative evaluations of a method for estimating recharge 
to the desert basins of Navado,. J.Hydrol., Vol.31, 335-357. 

Yakowitz, S. , and Duckstein, L. , 1980. Instability in 
aquifer identification: Theory and case studies, Uater Resour. 
Res., 16(6), 1045-1064. 

Yeh, U.U.G., 1970. Non-steady flow to surface reservoir. 
Journal of Hydr . Div. ASCE , 96(HY3), Proc. Paper 7137, 609-618. 

Yeh, U.U.G. , 1975. Aquifer parameter estimation. Journal 
of Hydr. Div. ASCE 101(HY9), 1197-1209. 

Yeh, U.U.G. and Tauxe, G.U., 1971. A proposed 
technique for unconflned aquifer parameters, J.Hydrol., 
Vol.12, No. 2., 117-128. 

Yeh, U.U.G. , 1986. Review of parameter identification 
procedures in groundwater hydrology. The Inverse problem, 
Uater Resour. Res., 22(2), 95-108. 

Yeh. U.U.G. and Y.S. Yoon, 1981. Parameter 
identification with optimum dimension in parameterization, 
Uater Resour. Res., 17(3), 664-672. 

Yeh, U.U.G., Y.S. Yoon, and K.S. Lee, 1983. Aquifer 
parameter Identification with krlging and optimum 
parameterization, Uater Resour. Res., 19(1), 225-233. 

Young, R. A. and Bredehoeft, J.D., 1972. Digital computer 
simulation for solving management problems of conjunctive 
groundwater and surface water systems, Uater Resour. 
Res., 8(3), 533-556. 

Yu, U. and Haimes, Y.Y., 1974. Multiple Optimisation for 
conjunctive use of ground and surface water, Uater Resour. 
Res., 10(4), 625-636. 



APPENDIX A 


An algorithm has been developed to process the water table 
head values by least square polynomial approximation. The salient 
features of this algorithm to calculate the least square 
polynomials are presented below: 


1. It normalizes the coordinates prior to the calculation of 
the coefficients of least square polynomials. 


2. It detects and deletes the terms of the polynomials which 
do not contribute significantly towards expanding power of 
the polynomial. This is done on the basis of students ’ t* 
test. 


3. 

4. 


5. 


6 . 


7 . 


It detects and deletes outliers based on the standard 
residue criterion. 


After deleting the terms and the outliers of 
square polynomial is re-computed with reduced 
terms and/or with reduced number of data points. 
Is continued till the ' t ’ test 


process 

warrant any further detection 
standard residues are acceptable. 


the least 
number of 
This 
does not 


of the terms and all 


After deciding the polynomial of minimum degree with 
minimum number of terms, corresponding least square 
coefficients are computed and out. 

It also calculates, if desired mean water table elevation 
(by the process of integration) and mean boundary gradient 
(by the process of differentiation) corresponding to 
lateral flow and boundary recharge respectively. 

The form of the least square polynomial and the 
coefficients are stored for further calculation. 


The flowchart of the algorithm is presented in the next 


page . 
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APPENDIX B 


The algorithm described in Appendix A provides the 
coefficients of polynomials at each well point. These 
coefficients help in the computation of spatial and temporal 
derivatives of the water table head at the space point under 
consideration (i.e. the space points for which the parameters are 
estimated) in each of the time periods. The following steps are 
involved. 


1 - It reads the coefficients of the least square polynomials. 

2. It reads the rainfall recharge values at all the 

observation points and for all the time periods. 

3. It reads the groundwater withdrawals, evaporation and 

canal seepage values for all space points and time 
periods . 

4. It reads the initial values, constants, coordinates of 

the observation points and its area. 

5. It calculates the second spatial derivatives and first 

temporal derivatives of head at all observation points. 

6. It assigns the recharge and discharge values for the 
corresponding observation points. 

7. An objective function is called through a sub-program and 
constraints and initial values are calculated. 

8. The objective function inturn calls a nonlinear 
optimization routine (SUMT) to estimate optimal estimates 
of aquifer parameters. 

A flow chart is presented In the next page indicating the 

algorithm used for the aquifer parameter estimation. 




FLOW CHART FOR AQUIFER PARAMETER ESTIMATION 






APPENDIX C 


Estimation of water table heads from a given set of data, 
such as transmissivity, storage factor, withdrawal rates etc., is 
called as direct problem in groundwater hydrology. Prickett and 
Lonnquist (1971) developed an algorithm for prediction of water 
table heads using finite difference method. In their algorithm 
they incorporated Iterative Alternating Direction Implicit method 
(Peaceman and Rachford, 1955). It involves first, for a given 
time increment, reducing a large set of simultaneous equations 
down to a number of small sets. This is done by solving the node 
equations by Gauss elimination by keeping individual columns/rows 
constant. The set of column/row equation is then implicit in the 
direction along the column/row and explicit in the direction 
orthogonal to the column/row alignment. 

The process of calculating heads along columns or rows 
in the digital model includes first computing the values of G and 
B as explained by Prickett and Lonnquist (1971U for the nodes of 

a column or row in increasing order. In the course of 

; 

accomplishing this, the head at the last node of the column or 
row is found. After completing the calculation of heads in an 
individual column or row, the computer proceeds to the next 
column or row until all have been processed satisfactorily. 



In the present study, after estimating transmissivity and 
storativity by inverse problem the same are cross checked by 
direct problem. Direct problem has also been used to predict 
water table heads after introducing a hypothetical drainage 
channel in the study area. 

A 20 X 19 grid was considered for the simulation of 
heads. The input data consists of heads at start of time 
increment, storage factor, withdrawals, transmissivity, number of 
time increments etc.. Storage factor (SFj) is calculated as 

SFj = 7.48 * S * Ax * Ay CA.l) 

where S = storativity, 

dx,A y = grid interval in x, y directions 


The error is estimated as 

Error = Q * DELTA / 10 * SF1 (A. 2) 

where Q = withdrawals, 

DELTA = time/ S (A. 3) 

where 6 is a function of time step. 


The process associated with this algorithm is 
presented In the next page. The program can do the following 

operations : 

1. It reads parameter and default value cards for all the 
space points. 

After reading node cards it. makes prediction of heads. 


2 . 




ADVANCE TIME, UPDATE H^ , AND MAKE PREDECTION 


SET SUM OF CHANGES IN HEAD, TO ZERO AND 
ADVAN CE I TERATION 


COMPUTE B AND G ARRAYS FOR A COLUMN 



COMPUTE HEADS FOR COLUMN AND SUM THE CHANGES 
IN HEAD INTO E 


ALL COLUMNS PROCESSED ? 


COMPUTE B AND G ARRAYS FOR A ROW 


COMPUTE HEADS FOR ROW AND SUM THE CHARGES 
IN HEAD ' INTO E 


ALL ROWS PROCESSED ? 


'CONVERGENCE (E ERROR) 1 


PRINT RESULTS 



FLOW CHART FOR DIRECT PROBLEM (AFTER 
PRICKETT AND LONNQUIST 1971 ) 










3. Uith advance in time, it computes B and G arrays for a 
column and make changes in heads. 

4. After processing all the columns, B and G arrays are 
estimated with respect to a row and changes are made in heads. 

5. Uhen heads are re-estimated for all the rows also it 
checks for convergence. 

6 . Uhen the convergence is satisfied it prints predicted 
head values for this time increment. 

The above steps are repeated till 
int ervals . 


7. 


the end of time 



